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Summary. We present a two-dimensional (2D) mathematical model of a highly 
concentrated suspension or a thin film of the rigid inclusions in an incompressible 
Newtonian fluid. Our objectives are two-fold: (i) to obtain all singular terms in 
the asymptotics of the overall viscous dissipation rate as the interparticle distance 
parameter 8 tends to zero, (ii) to obtain a qualitative description of a microflow 
between neighboring inclusions in the suspension. 

Due to reduced analytical and computational complexity, 2D models are often 
used for a description of 3D suspensions. Our analysis provides the limits of validity 
of 2D models for 3D problems and highlights novel features of 2D physical problems 
(e.g. thin films). It also shows that the Poiseuille type microflow contributes to a 
singular behavior of the dissipation rate. We present examples in which this flow 
results in anomalous rate of blow up of the dissipation rate in 2D. We show that 
this anomalous blow up has no analog in 3D. 

While previously developed techniques allowed to derive and justify the leading 
singular term only for special symmetric boundary conditions, a fictitious fluid ap- 
proach, developed in this paper, captures all singular terms in the asymptotics of 
the dissipation rate for generic boundary conditions. This approach seems to be an 
appropriate tool for rigorous analysis of 3D models of suspensions as well as various 
other models of highly packed composites. 

Key words: concentrated suspensions, overall or effective viscous dissipa- 
tion rate, Stokes flow, discrete network approximation, variational bounds, 
Poiseuille flow. 



1 Introduction 



Many classical and novel engineering processes involving multiphase flow re- 
quire to capture the overall behavior of suspensions. The problem of the be- 
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havior of suspensions is important in geophysics (mud-flow and debris flow 
rhcology), pharmacology (drugs design), ceramics processing among others. 
Wide range of experimental (see, e.g., [1,3,4,20,31,37,44-46,52]), and nu- 
merical (see, e.g., [17,35,40,49,53,55]) results are available. 

Two-dimensional (2D) models of three-dimensional (3D) suspensions are 
often used in numerical simulations (e.g. [21, 34, 47]), because they require less 
computational effort. In addition, 2D suspensions can describe biological thin 
films [56] and certain types of 3D concentrated suspensions of uniaxial thin 
rods (see e.g. [39]). 

We consider a 2D mathematical model of a non-colloidal (Browning mo- 
tion can be neglected) concentrated suspension of neutrally buoyant rigid 
particles (inclusions) in a Newtonian fluid. The suspension occupies a 2D do- 
main fi, and rigid inclusions are modeled by disks of equal radii. The main 
objective is to characterize in a rigorous mathematical framework the depen- 
dence of the effective rheological properties (e.g. effective viscosity, effective 
permeability, effective viscous dissipation rate) of such suspension on the ge- 
ometry of inclusions array and applied boundary conditions on df2. We focus 
on highly packed suspensions when the concentration of inclusions is close to 
maximal, which means that the distance between neighboring inclusions (in- 
terparticle distance) is much smaller than their sizes. We consider an irregular 
(non-periodic) array of disks and our analysis takes into account the variable 
distances between adjacent inclusions. 

The main features of the problem under consideration are the high concen- 
tration of the inclusions and the irregular geometry of their spatial distribu- 
tion. The key quantity of interest in describing the effective rheologicahprop- 
erties of suspensions is the rate of viscous dissipation of energy W — W(u), 
where u is the velocity of the incompressible fluid (see the precise definition 
in Section 2). 

Initially our interest was motivated by the problem of dewatering process 
which led to the sedimentation due to gravity problem [19] arising in the dis- 
posal of solid waste in industrial and municipal pollution management. In this 
problem particles are close to touching. In particular, shear external boundary 
condition was brought up by experimental and numerical studies [27,29,30] 
and it was addressed mathematically by analysis of 2D model in [10]. This 
led us to the question of limits of validity of 2D models for description of 3D 
suspension. 

Note that the problem of sedimentation is quite difficult in mathematical 
content and far from being settled. Even in the case of the dilute limit this 
problem is not completely understood mathematically and we mention here 
the number of recent mathematical studies in this direction [7, 51]. 

For highly concentrated suspensions of rigid inclusions, W exhibits a sin- 
gular behavior (see e.g. [8,10,24,28,50,53]) and its understanding is a fun- 
damental issue. A formal asymptotic analysis of the singular behavior of the 
viscous dissipation rate in a thin gap between a single pair of two closely 
spaced spherical inclusions in a Newtonian fluid was performed in [24] . In this 
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work only translational motions of inclusions but not rotational were consid- 
ered and the asymptotics of the form CS~ X +0(ln 1/8), where 5 is the distance 
between two spheres, was obtained. Based on the analysis of a single pair of 
spheres the authors of [24] suggested that the asymptotics of the effective vis- 
cosity of the periodic array of inclusions is of the same form, that is, the main 
singular behavior is of 0(5^ 1 ). Next, a periodic array of inclusions in a New- 
tonian fluid was considered in [50]. Under the assumption that all inclusions 
follow the shear motion of the fluid (formula (5) in [50]) it was shown that 
W = 0(S" 1 ). This assumption is analogous to the well-known Cauchy-Born 
hypothesis in solid state physics, which is known to be not always true [25]. 
Indeed, in the case of suspensions it was shown in [10] that for shear external 
boundary conditions the inclusions may not follow the shear motion. More- 
over, it was shown in numerical studies of [53] that the asymptotics 0(S^ 1 ) 
may or may not hold for suspensions of a large number of inclusions with 
generic boundary conditions. There it was observed numerically that while in 
some cases the asymptotics of the effective viscosity is of order 1/(5, in other 
cases it is of order In 1 /5. Also the problem of the exact analytical form of the 
singular behavior for generic suspensions was posed in [53] (p. 140) which mo- 
tivated subsequent studies of [8, 10] and the present paper, where a network 
approximation approach was developed for this type of problems. 

In [8] a formal asymptotic analysis of the effective viscosity in 3D for 
a disordered array of inclusions was performed. In a view of discrepancies 
between predictions of the formal asymptotic analysis [24] and numerics [53] 
such formal asymptotics requires a mathematical justification. In [8] for special 
(extensional) boundary conditions the leading term of the asymptotics of the 
effective viscosity as 5 — > 0, where 5 is the characteristic spacing between 
neighboring inclusions, was justified in a 2D model. In subsequent work [10] 
this leading term, that exhibits a so-called strong blow up of order (5~ 3 / 2 , 
was analyzed. It turned out that in many important cases, e.g. when shear 
boundary conditions are applied, this term degenerates, so the next term of 
order <5 -1 / 2 , that exhibits a so-called weak blow up, becomes the leading term 
of the asymptotics in many physical situations. 

However, the techniques of [8] are only capable to capture the strong blow 
up but not the other singular terms in the asymptotic expansion of the ef- 
fective viscosity. By contrast, a fictitious fluid approach introduced in this 
paper allows to derive the complete asymptotic expansion of the overall vis- 
cous dissipation rate in which all singular terms are captured and justified. In 
particular, we ruled out singular terms other than presented in Theorem 2.1. 
Previously [8-10,24,28,50,53] inclusions in a dense suspension were charac- 
terized by the sets of their translational and rotational velocities, so-called 
discrete variables. Our analysis shows that in order to obtain the complete 
asymptotics of singular behavior it is necessary use an additional set of discrete 
variables, permeation constants. To the best of our knowledge these discrete 
variables have not been used in previous studies of dense suspensions. We now 
explain the physical consequences of this asymptotics. 
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The key feature of rheology of concentrated suspensions is that the dom- 
inant contribution to the overall or effective viscous dissipation rate (and 
therefore, to the effective viscosity) comes from thin gaps (lubrication re- 
gions) between closely spaced neighboring inclusions [43]. The mathematical 
techniques introduced in the above mentioned works [8, 24, 50] utilized this 
observation. More specifically, they took into account certain types of relative 
movements of inclusions which resulted in the corresponding microflows in 
the gaps. The formal asymptotics in [24] was based on analysis of the squeeze 
motion, when two inclusions move toward each other along the line joining 
their centers (see Fig. 2.7c) but did not provide the detailed analysis of other 
motions, which was sufficient for certain type of boundary conditions (e.g. ex- 
tensional boundary conditions) but not sufficient for others. A justification of 
the formal asymptotics was not considered in [24] (see also [28] where similar 
results were obtained). 

In [8] four types of relative motions of neighboring inclusions were consid- 
ered which result in a singular behavior of dissipation rate: the squeeze (Fig. 
2.7c), the shear (Fig. 2.7b) and two rotations (Fig. 2.8). While it was sufficient 
for the analysis of the leading term of the overall viscous dissipation rate (in a 
suspension of free particles) which was the goal of [8] , in the present paper we 
observed that this analysis does not provide a complete picture of microflows. 
Indeed, the Poiseuille flow in 2D also contributes into the singular behavior. 
Examples in the present paper suggest that when an external field is applied 
to inclusions in a suspension, this Poiseuille flow may give rise to an anoma- 
lous strong rate of blow up (called a superstrong blow up, of order 5~ 5 ^ 2 ) of 
the viscous dissipation rate, whereas for suspensions of free inclusions there 
is at most strong blow up (of order <5~ 3 / 2 ). 

As explained above the previous studies [8, 24, 28, 50] provided only a par- 
tial analysis of the single behavior of the viscous dissipation rate and therefore 
did not provide the complete physical picture of microflows. In this work we 
obtain the complete asymptotic description of the singular behavior of viscous 
dissipation rate which led us to the complete description of microflows (Fig. 
2.7-2.10). 

The techniques of [8, 10] and the present paper are based on the discrete 
network approximation. Discrete networks have been used as analogues of the 
continuum problems in various areas of physics and engineering for a long 
time (see, for example, [12] and references therein, and see also the recent re- 
view [48] for various applications of networks in social and biological studies). 
However, the fundamental issue of relationship between a continuum problem 
and the corresponding discrete network was not addressed until recently. 

The pioneering study of electro-magnetic properties of high contrast 
materials arising in imaging using a rigorous network approximation was 
done in[14-16]. There the local resistivity of a periodic medium with con- 
tinuously distributed properties was given by the Kozlov's function ([38]): 
p(x) = Const e~ s ( x ^ e , where S(x) is a periodic smooth phase function, e is 
an aspect ratio of the material properties, a small parameter of the problem 
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(as e — > 0, p(x) describes a resistivity of a high contrast composite). The effec- 
tive resistivity obtained and justified there was given in terms of the principal 
curvatures k + , kT of the function S(x) at the saddle point of S(x). 

The discrete network approximation for a medium with piecewise constant 
characteristics, which correspond to particle-filled (particulate) composites, 
was developed in [12]. In this case, the principle curvatures and saddle points 
of the phase function are not determined. In [12] an infinite contrast material 
(e = 0) with irregularly or non-periodically distributed circular fibers was 
considered and asymptotics of the effective conductivity of such a composite 
as 5 — > was derived and justified. In [12] the Keller's asymptotic solution 
[36] for two closely spaced disks was used as a building block for construction 
of the network for a large number of inclusions (disks). The network model 
is a linear system on a graph with appropriate boundary conditions. The 
sites or vertices of the graph correspond to the centers of inclusions {xi, i = 
1, . . . ,7V} and the values of the local fluxes between adjacent inclusions Xi 

and Xj described by the Keller's formula: n^J^, where R is the radius of 

disks, Sij is the distance between the closely spaced neighbors, are assigned 
to the edges of the graph. The notion of the neighbors is introduced via the 
Voronoi tessellation. The asymptotic closeness of the effective conductivity of 
continuum and discrete problems was justified by using the direct and dual 
variational bounds matching to the leading order which was obtained to be 
of order as 5 — > 0, where 5 = max Sij. 

In [13] the notion of the characteristic interparticle distance was general- 
ized for a broad class of 2D geometrical patterns (highly non-uniform arrays 
of inclusions). For this class of arrays an explicit error estimate for the dis- 
crete network approximation was obtained. A scalar conductivity problem 
analogous to [12, 13] in 3D was considered in [11] where both asymptotics of 
the effective conductivity and its error estimated were obtained. While the 
formulation of the problem is analogous to the 2D one considered in [12, 13], 
the connectivity patterns of conducting inclusions are much more complex 
in 3D. That is why the justification of this approximation and derivation of 
its error estimate required introduction of new techniques (e.g., the central 
projection partition). By contrast, our work shows that 2D vectorial problems 
may exhibit features which have no analogs in similar 3D problems. 

Next the method of the discrete network approximation was extended to 
a vectorial problem [8] , where formal asymptotics for 3D suspensions was de- 
rived. The approach developed in [12, 13] for the justification of the formal 
asymptotics in the scalar problem could not be readily applied to the vec- 
torial problems, which exhibit new features and require new technical tools. 
In particular, the key point of discrete network approach is a construction of 
trial functions in the variational formulation of the problem for determining 
the effective properties of composites. In vectorial problems (unlike in scalar 
ones) such trial functions (for the direct variational principle) must satisfy 
the divergence free condition. Also the dual variational formulation involves 
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integral constraints - balances of linear and angular momenta. The balance 
equations on an inclusion involve all its neighbors that have their own neigh- 
bors etc. This accounts for a long-range hydrodynamic interactions and shows 
that the pairwise (local) analysis (such as [24, 28]) is not sufficient and a global 
analysis is necessary for the whole array of inclusions. That is why a general- 
ization of scalar techniques from [12, 13] was done for the leading term only 
in the asymptotics of the effective viscosity in the 2D model [8] . It seems that 
a straightforward generalization of the approach in [8] becomes increasingly 
difficult and probably impossible due to the global nature of above mentioned 
constraints. 

The objectives of this work are two-fold: (i) develop a new approach, 
called a fictitious fluid approach, that allows to effectively deal with global 
constraints (balance equations and the divergence free condition), (ii) apply 
the fictitious fluid approach to derive and justify all singular terms in the 
asymptotics of the overall viscous dissipation rate. Such an asymptotic for- 
mula results in a complete description of microflows in suspensions, as oppose 
to the partial description of the microflow in [8, 24, 28, 50]. 

The techniques in [11-13] for scalar problems in both 2D and 3D were 
based on a direct (one-step) discretization of the original continuum problem. 
In vectorial problems [8] the direct discretization allowed to obtain only the 
leading term for the extensional viscosity (strong blow up) but not for the 
shear viscosity (weak blow up) due to the global constraints. In the present 
work we use a two-step discretization procedure based on the fictitious fluid 
approach to derive and justify both strong and weak blow up terms. 

We briefly describe now the idea behind the fictitious fluid approach. As 
mentioned before, it consists of two steps. In step 1 we introduce a fictitious 
fluid domain which comprises thin gaps between neighboring inclusions. The 
dissipation rate restricted to this domain is denoted by Wn- We show that 
for generic Dirichlet boundary conditions it describes the singular behavior in 
the following sense: 

W = W n + h.o.t. as 6^0. 

In steg 2 we perform a discretization of Wn, that is, the continuum prob- 
lem for Wn is reduced to an algebraic problem on a graph, called a network 
problem. By using the fictitious fluid approach in this step most of the diffi- 
culties are eliminated. The network problem is a minimization of a quadratic 
form whose coefficients depend on R, fj, and boundary data. The quadratic 
form on the minimizing set of discrete variables is called the discrete viscous 
dissipation rate and denoted by X. Unknowns of this problem are vectors 
U = {U l }f =11 ijj = {to l }f =1 , the translational and angular velocities of inclu- 
sions, respectively, and a collection of numbers (3 — {/%} characterizing the 
Poiseuille microflow between a pair of inclusions. The discrete dissipation rate 
X is given by: 

J = Ti(/3)(T 5 / 2 +Z 2 (U,u;,/3)<r 3/2 + Z 3 (U, ^(T 1 / 2 , as 5 -> (1.1) 
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where k = 1,2,3, are explicitly computable quadratic polynomials of U, 

As a result we obtain the following asymptotic formula for the generic 
Dirichlet boundary conditions: 

In fact, we prove the following result about the error term: 

\W - 1\ < n (^p E^fcj) + E 2 (C/ 4 - U j ) + E 3 (u/) j 

where Ei, E2, E3 are quadratic polynomials of fy, difference U l — U 3 and u> 1 , 
respectively, whose coefficients are independent of 6. 

Finally, we analyze the physical ramification of the obtained asymptotic 
formula by presenting several examples. We construct an a example of a sus- 
pension in a strong "pinning" external field, where W is of order <5~ 5 / 2 (super- 
strong blow up). We also show an example of the superstrong blow up due to 
the boundary layer effect. Note that to the best of our knowledge this rate of 
blow up was not observed before and we call it an anomalous rate. For generic 
suspensions (free particles or a weak external field) we expect that (3 = 0(5) 
and, therefore, the first and the third terms of (1.1) are of the same order 
0((5 -1 / 2 ). For a hexagonal array of inclusions we prove that W exhibits the 
strong blow up (of order <5~ 3//2 ) and (3 = 0. Note that a typical close packing 
array in 2D is "approximately" hexagonal. 

The paper is organized as follows. In Section 2 we give a mathematical 
formulation of the problem (Subsection 1), describe the fictitious fluid ap- 
proach, and present our main results in Theorems 2.1, 2.2 (Subsection 2). In 
Subsection 3 we construct our discrete network and discuss how local flows in 
thin gaps between neighbors (microflows) contribute the overall viscous dissi- 
pation rate and state the theorem about a representation of the error term of 
the discrete approximation. In Section 3 we discuss main and present exam- 
ples. Section 4 is devoted to the the fictitious fluid problem. In Section 5 we 
present results related to our discrete network. In section 6 prove two lemmas 
about approximation of the overall dissipation rate by the dissipation rate of 
the fictitious fluid and about approximation of the fictitious fluid (continuum) 
dissipation rate by the discrete one. Conclusions are presented in Section 7. 
The proofs of some auxiliary facts are given in Appendices. 
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2 Formulation of the Problem and Main Results 
1 Mathematical Formulation of the Problem 

Consider an irregular or non-periodic array of N identical circular disks B l , of 

N 

the radius R distributed in a rectangular domain Q. Denote by Qf = [J B % 

i=i 

the fluid domain which is occupied by incompressible fluid with viscosity \x 
(see Fig. 2.1). Disks B % represent absolutely rigid inclusions. Inertia of both 
fluid and inclusions is neglected. In the fluid domain flp consider the following 
boundary value problem: 

(a) (j,Au — Vp, x £ ftp 

(b) V -11 = 0, x£fl F 

(c) u = U l + Ruj\n\e 2 - u|ei), x £ dB\ i = l...N 

(d) / cr(u)n l ds = i = l...N 

Job* 

(e) / n' x er(it)n 4 ds = 0, 

J dB* 



(2.1) 



i = l...N 

x£dn 

where u(x) is the velocity field at a point x £ Op, p(x) is the pres- 



dB 

[(/)« = /, 




Fig. 2.1. Domain Of occupied by the fluid of viscosity fi, and disordered array of 
closely spaced inclusions B l 



1 / du ' On ■ 

sure, eriu) = 2llD(u) — pi is the stress tensor, Dij(u) = - — — - + — — 

2 \ oxj OXi j 
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i, j = 1, 2, is the rate of strain which satisfies the incompressibility condition: 
tr D(u) — 0, another form of (2.16). The vector n 1 = (n\,ni£) is the outer 
normal to B % . Constant vectors U l — {U\, f7|) and scalars uj 1 , i — 1, . . . ,N, 
which are translational and angular velocities of the inclusion B l , respectively, 
are to be found in the course of solving the problem. 

Here N is closed to maximal packing number N max = N max (ft, R). This 
number is finite and \N — N max \ depends on the small parameter S called 
interpaticle distance rigorously defined below in Definition 1. 

We consider the linear external boundary conditions of the form: 

'-*-(:-.)(;) (2 - 2) 

where the components a, b, c of the matrix A are given constants. Note that 
the most general form of the linear boundary conditions is / = / + Ax where 
/ is a constant vector. Observe that when a = and b = — 1/c the vector 
Ax corresponds to a rotation and f to a translation of the boundary, hence, 
/ describes the rigid body motion of dfl. Hereafter, we exclude this trivial 
motion from our consideration assuming that a — and b = —1/c does not 
hold simultaneously in (2.2) and / = 0. We use such boundary conditions 
for two reasons: a) for technical simplicity, which does not lead to the loss of 
generality; b) they include the shear (when a = c = 0, b = 1) and extensional 
(a = 1, b = c = 0) boundary conditions (see e.g. [10]) which model two basic 
types of viscometric measurements (see Fig. 2.2). It is possible to extend our 
results to arbitrary Dirichlet boundary conditions / S H 1 ^ 2 (df2) satisfying 

/ / • nds = 0. 

Jan 



J—L 1 LA 




\ \ t / / 

(b) extensional 



(a) shear 

Fig. 2.2. Shear and extensional external boundary conditions 



For an arbitrary set A C f2 F consider the following integral: 



W A (v) = - I o-(v) : D(v) = n I D(v) : D(v)dx 

2 J A _ J A 



/' 



dvi 
dx 



A 

dvi 
dy 



dv 2 
dx 



dv 2 
dy 



dx 



(2.3) 
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where v = (v\, v 2 ). Then the variational formulation of (2.1) is: 

Find u e V, such that Wq f {u) = minWo F (v), (2.4) 

i)6 V 

where the set of admissible vector fields V is defined by 

V = {v e H 1 (f2 F ) : V • v = in Q F , v = f on dC2, 

1 (2.5) 
v = U l +u> i x(x- x r ). x e dB\ i=l,...,N}, 

Wq f (u) is called the (continuum) viscous dissipation rate [42] and it is the 
principal quantity of interest in the study of overall properties of suspensions. 
We will use the following notation: 

W:=W Qf {u). (2.6) 

The key feature of our problem is that we study suspensions where concentra- 
tion of inclusions is close to its maximum. Therefore, the domain Q F depends 
on the characteristic intcrparticle distance parameter 5. Our main objective 
is to derive and justify an asymptotics of W as S — > 0. We will show that 
the coefficients of this asymptotic formula are determined by the solution 
to a discrete network problem, which determine the discrete viscous dissipa- 
tion rate. As explained in Introduction, W determines a number of important 
measurable properties of suspensions. 



2 The Fictitious Fluid Approach and Discretization 

In Appendix 1 we show that using the notion of Voronoi tessellation we can 
decompose the domain ftp into necks II and triangles A: Q F = II LI A (see 
Fig. 2.4). Necks connect either two disks (Fig. 2.3a) or a disk and a part of 
the boundary dfl called a quasidisk (see Fig. 2.36), that is, necks connect 
neighbors. The velocities of quasidisks are given by the prescribed boundary 
conditions (2.2). Note that near the boundary when quasidisks are involved 
the "triangles" are actually trapezoids. With slight abuse of terminology, we 
also call them triangles. 

We distinguish boundary disks (quasidisks) and interior disks and intro- 
duce two sets of the corresponding indices. For indices of interior disks we use 
the notation I = {1, . . . , N}. If the disk B l centered at Xi is a quasidisk then 
i belongs to the set B of the indices of quasidisks. Also denote by A/i the set 
of indices of all neighbors of B\ 

For a given array of the disks and quasidisks B % centered at x i} the dis- 
crete network is the graph Q = (X, £), with set of vertices X = {xi : i e 
IUB} and set of edges £ — {eij : i £ I, j £ Ni} with connecting neighbors 
B i and BK 

As mentioned in Introduction our main approach in study of the asymp- 
totics of W as S — > consists of two steps. This two-step approach allows to 
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Fig. 2.3. (a) neck connecting two disks; (6) neck connecting disk B l and quasidisk 



separate geometric construction of the network and its subsequent asymptotic 
analysis. 

In the first step, we show that the minimization problem (2.4) in the fluid 
domain ftp can be approximated by a "fictitious fluid" problem in which fluid 
is assumed to occupy necks U = [J llij between closely spaced neighbor- 

ing inclusions (the shadowed region in Fig. 2.4). We call 77 the fictitious fluid 
domain. On the boundary of the complementary part of the domain (triangles 
in Fig. 2.4) the relaxed incompressibility conditions: 



v ■ nds — 0, i € I, j, k € A/i, 



(2.7) 



are imposed. This reflects a well-known physical fact that for densely packed 
suspensions the dominant contribution to the viscous dissipation rate over the 
fluid domain comes from those necks. 

Below we show that the functional W is decomposed as follows: 



W = Wu + Wa, 



(2- 



where Wn is the overall viscous dissipation rate of the fictitious fluid defined 
below in (2.10) and Wa is the remainding contribution from the domain A. 

Consider the problem of minimization of the functional Wn, defined by 
(2.3) over the fictitious fluid domain II, in the following class of functions: 

V n = {v € H x {n) : V- v = in II, / v ■ nds = for all A ijk e A, 

Jd/\ iJk 

v = U l + Ruj\n\e 2 - n\ex), x £ dB\ i = l,...,N, v = f on dll n df2} , 

(2.9) 
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•i r 


V 7 ) 
~~J k J 













Fig. 2.4. The decomposition of the original fluid domain Of into the fictitious fluid 
domain and the set of triangles 

where U l and lo 1 , i — 1, . . . , N, are arbitrary constant vectors and arbitrary 
constants, respectively. The zero-flux condition through dAy^ (2.7) is inher- 
ited from the problem in the original fluid domain fip due to incompressibility 
condition of the fluid in triangles. Such a condition is a necessary (but not 
sufficient) for V • v — in the triangle Ay^. 

We define the overall dissipation rate of the fictitious fluid by 

W n = min W n {v). (2.10) 

veV n 

The first principal result of this paper is that the dissipation rate W can 
be approximated by the rate Wn and Wa can be neglected. To show this we 
need to introduce a small parameter of the problem, which is a characteristic 
interparticle distance 5. 

For each pair of neighbors x%, Xj define 

{la;, —Xj\ — 2R, when i, j e I, 
\xi — Xj \ ~ R, when either i or j G B, (2-11) 
\xi —Xj\, when i, j e B. 

As mentioned above, we study domains with closely spaced neighboring disks. 
More precisely, for all pairs of neighbors we assume that following close- 
packing condition holds. 

Definition 1. Write the minimal distance 8ij (see Fig. 2.3) between any two 
neighboring disks B % and B° in the form Sij = Sdij, where dij is such that 
< c\ < dij < C2 for some absolute constants c\, C2- // the parameter 5, called 
the characteristic interparticle distance, tends to zero, then Qp is said 
to satisfy the close-packing condition. 

We remark that this definition describes uniformly dense arrays of disks. 
A more general definition which covers a notion of a "hole" corresponding to 
the void space in the composite is introduced and discussed e.g. in [11, 13]. 
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Hereafter we call the array of inclusions under consideration a quasi- 
hexagonal array (e.g. in [13] it is referred to as "randomized hexagonal"). 
Recall, that for such arrays all neighbors are closely spaced and a typical 
number of nearest neighbors for a disk is six. 

Hence, the mathematical thrust of step one of the fictitious fluid approach 
is in showing that the overall rate of the energy dissipation Wn of the fictitious 
fluid captures the singular behavior of W, defined by (2.6), as S — > 0. More 
precisely, in Section 5 we prove the following theorem: 

Theorem 2.1 (Approximation by the fictitious fluid). Suppose an array 
of inclusions satisfies the close packing condition. Let W be the overall viscous 
dissipation rate defined by (2.5)-(2.6) and Wn be the viscous dissipation rate 
of the fictitious fluid defined by (2.10). Then the following asymptotic formula 
holds: _ _ 

\ W - Wn \ = Oi?l*) as 5^0. (2.12) 



In step two, we study asymptotics (blow up) of the overall viscous dis- 
sipation rate W asj5 — > 0. In view of step one, this is reduced to finding 
of asymptotics of Wn- The latter is done by a construction of a discrete 
network approximation and introduction of a so-called discrete viscous dissi- 
pation rate X. To show closeness of the continuum and the discrete dissipation 
rates, Wn and X, respectively, we employ the direct and dual variational tech- 
niques [8,12,13]. 

Also note that the conditions (2.1d,e) in the original problem led to sig- 
nificant technical difficulties in variational analysis of the overall viscous dis- 
sipation rate of [8], which is why the analysis of [8] is restricted to its leading 
singular term. In contrast, analogs of these conditions in the fictitious fluid 
problem are satisfied automatically by construction, which results in substan- 
tial simplification of the analysis and thus allows to capture all singular terms. 

The approximation of the overall viscous dissipation rate W by the discrete 
dissipation rate X is given by the following theorem. 

Theorem 2.2 (Approximation of the continuum dissipation rates by 
the discrete one). Suppose Qf satisfies the close packing condition and 
X= min Q(U, u>, /3), where the positive definite quadratic form 

(U,u>,j3)€"R. 

Q(U, w, 0) = Ii (/3)<T 5 / 2 + X 2 (U, w, /3)<T 3/2 + X 3 (U, u, 0) J" 1 / 2 

on the class of admissible discrete variables (U, u>,/3) € 1Z defined in (2.33)- 
(2.37). Then the following approximation to the viscous dissipation rate holds: 

\K^I1 = 0(5^) as 6^0, (2.13) 



Remark 2.1. Zj (i = 1, 2, 3) are explicitly computable quadratic polynomials of 
(U, w,/3) € TL from equations (2.36)-(2.37) below. These polynomials depend 
only on boundary data /, viscosity a and geometry of Q F . 



14 



Leonid Berlyand, Yuliya Gorb, and Alexei Novikov 



The next subsection is devoted to the construction of these quadratic poly- 
nomials. We first introduce a set of discrete variables (U, u>,/3), define a new 
variable 0ij in each neck 11^ and explain how the quadratic form Q is ob- 
tained. We will also discuss the underlying structure of the flow in a neck 
(microflow) and physical ramifications of the obtained asymptotics (2.13) and 
(2.38) below. 





3 Construction of the discrete network. Microflows. 

(a) We begin with the discretization of the boundary conditions. Denote U = 
{L^jigiuB, w = {a/}i e iuB where on the boundary <9I2, that is, for these 
velocities are given by boundary conditions as follows: 



Xi e dn lat 

IT={)\ and uj 1 = 7 1 " Q "^_ (2.14) 

Xi e <)<> : , h 

where Xi is the center of the quasidisk B l and dfl + and df2~ are the upper 
and lower parts of dfl, respectively, and dQ^ at and dfi^ at are the lateral (left 
and right) boundary. 

(b) Discretization of the incomprcssibility condition is implemented as fol- 
lows. Decompose the domain Q F into curvilinear hexagons Aijk as in Fig. 
2.5: 

= [J Aijk- 

Each Aijk consists of the line segments lij, £jk, Iki (Fig. 2.5a) connecting 
disks £>*, B- 7 ', B k and arcs a^, a^, of the corresponding disk. 3 Then the 
weak incompressibility condition (2.7) for the class Vn (2-9) becomes 

/ u ■ nds = 0, for any Aijk, i <= I> J, k G A/",. (2.15) 

In order to continue our analysis at this point we must introduce a new 
set of discrete variables. Here we define permeation constants: 

fij ! = 15 / u ■ nds > i G ! ' j 6 A/ ' 4 ' 

fii = nl f-nds,eijCdf2(i,jeM), 



3 In the case when the disk B l has two quasidisk neighbors B j and B k , that is, 
when B k is in the corner as in Fig. 2.56, then the domain Aijk is actually a 
curvilinear pentagon. By slight abuse of terminology we still call it a "curvilinear 
hexagon" . 
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(a) (b) 

Fig. 2.5. (a) Decomposition of Of into curvilinear hexagons Aijk and line lij 
connecting neighbors B % and B 3 , (b) Construction of Aijk at the corner of the 
boundary dO 

where is the line segment joining two neighbors B % and B 3 , u £ Vjj and 
n is an outer normal to Aijk- Then (2.15) can be rewritten as: 

Pii + P* 3 k +^i + ^J a U *- nlds + JiJ a W ■ nHs + JiJ a Uk ■ nkds = °> 

for i £ I, j, k £ JVi (n l is a unit outer normal to dB l ), which can be further 
simplified as 

/?*■ + (3* k + p* ki + {U l + U 3 )p 13 + (U 3 + U k )p> k + (U k + U' l )p kl = 0, (2.17) 

where vectors q v and p lJ are the unit vectors of the local system coordinate 
of two neighboring disks B l and B 3 as in Fig. 2.6. We call (2.17) the weak 
incompressibility condition. This formula explains the scaling in (2.16). In- 




Fig. 2.6. Unit vectors of the local coordinate system 
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deed, from dimensional analysis L u-nds must be divided by a characteristic 
lengthscale which is R in the curvilinear hexagonal Aijk- 

(c) We finally discretize the Stokes equations in necks. This discretiza- 
tion is based on lubrication theory [2,6,43]. This theory describes thin film 
flows between two solid bodies sliding relative to each other. It is well-known 
(see e.g. [43]) that in dense suspensions the dominant hydrodynamic contri- 
bution to the viscous dissipation rate occurs in necks between closely spaced 
inclusions, where lubrication equations are relevant. Lubrication theory de- 
termines fluid motion in such necks as a result of relative kinematic motion 
of the neighboring inclusions. This raises a question of the classification of 
microflows, that is, local flows in necks between two closely spaced neighbors. 

Recall three classical types of microflows between two parallel plates: the 
shear, the squeeze, and the Poiseuille flow. The last one, however, is not re- 
lated to motions of two plates relative one to another and therefore it is not 
described by lubrication theory. In this work we show that exactly these three 
types of microflows fully describe the motion of the fluid between two neigh- 
boring disks. Asymptotic analysis of microflows between two parallel plates 
technically is much simpler than that for inclusions with curvilinear bound- 
aries, which is needed for suspensions. 

As shown below, in classical 3D problem the Poiseuille microflow between 
two inclusions does not contribute to the singular behavior of the viscous dis- 
sipation rate. However, in analogous 2D problem (e.g. thin films) all three mi- 
croflows are present and, moreover, the Poiseuille flow results in anomalously 
strong singularity. Note that in the previous studies of the overall properties 
of suspensions the Poiseuille type microflow was not taken into account. 

Thus, it is necessary to analyze kinematics of a pair of neighboring inclu- 
sions when one moves relative to the other. To this end for a pair of neighbors 
B l and B- 7 ', centered at Xi and Xj, we choose the local coordinate system 
where the origin is at (xi + Xj)/2 and the y-axis is directed along the vector 
connecting Xi and Xj . 

For clarity of presentation we consider two interior disks only (that is, i e I, 
j € Mi fl I). An analogous construction in boundary necks [i e I, j e A/i fl B) 
is given in Appendix 2. 

There are exactly five elementary kinematic motions of inclusions. To see 
this we consider boundary conditions on dB l , dB J in (2.9) and first assume 
that u l — uji — 0. Then the conditions: u = U l on dB l and u — U J on dfP 
can be rewritten as follows: 

u\ dBi = \{IP + W) + \{U{ U() ei + \{Ui Ui)e 2 , 
and u\ dB3 =\{U l + U j ) - \{U{ - U{)e x - \{U\ - U{)e 2 . 

Hence, the translational velocities of disks are decomposed into three motions. 
First, when both inclusions and fluid move with the same velocity (see Fig. 
2.7a), and therefore this motion does not contribute to the singular behavior 



Anomalous blow-up of the Viscous Dissipation Rate 



17 




(a) parallel translation 




U'-ui) 



2(U/-U/) 

(b) shear 




(c) squeeze 



Fig. 2.7. Decomposition of translational velocities U l and U J into three elementary 
motions 



of the viscous dissipation rate. Second, when one inclusion moves relative to 
the other producing the shear type motion (Fig. 2.7b) and, finally, the squeeze 
type motion (Fig. 2.7c) of the fluid. 




§ (oj'+oj^inj e 2 -ri 2 e 1 ) 
§(a)'+aj j )(nje 2 +n 2 e l ) 




§(a) i -a>>)(n i 1 e 2 -ri 2 e 1 ) 
§(cx) i -o) J )(nje 2 +n i 2 e 1 ) 



(a) rotation in the 
same direction 



(6) rotation in the 
opposite directions 



Fig. 2.8. Decomposition of the angular velocities uj 1 and uj 3 into two elementary 
motions 



Similarly, assuming U l = U 3 = in conditions on dB l : dB 3 in (2.9) we 
decompose them into two relative elementary motions: rotations of the disks 
in the same direction (Fig. 2.8a), and rotations of the disks in the opposite 
directions (Fig. 2.86). Thus, 

R R 

u \dB, = — (uj 1 + uJ 3 )(n\e 2 - raid) + —(uj 1 - uj j )(n\e 2 ~ n\e{) 

R R 

and u\ dB j = — + tJ)(n\e2 +n|ei) + — (to 1 - aj j )(-n\e2 — n^ei), 
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where ri 1 = (n\, n|) is the outer normal to dB % . 

Next consider microflows corresponding to the four kinematic motions con- 
tributing to the singular behavior of the dissipation rate described above. We 
further choose the microflow ui u b which minimizes the viscous dissipation rate 
in the neck 77^ by imposing the natural boundary conditions on the lateral 
boundaries dll^j (Fig. 2.9). That is, the function ui u b minimizes Wn i3 in the 
class: 

Vij = {v e i7 1 (7Ty) : v • v = in 77^, v = g i on dB l , v = g } on dB j } 

(2.18) 

with g l = U l + Rw l {n\e 2 - n\e-y), g = U 3 + Rtu J (n{e 2 - n|e x ). (2.19) 



an; 



il 



Fig. 2.9. Boundary of the neck IJi 



Due to linearity the minimizer ui u b is decomposed into five vector fields 
corresponding to the relative motions of inclusions as in Fig. 2.7-2.8 as follows: 

uiub = \(U l + U*) + {(U l -IP). p«] ttl + [(U* -W)- g«]«a 
+ 73(u/ + LJ j )u 3 + Rfa* - uj j )u 4 , 

where functions u/-, k = 1,...,4 are minimizers of Wn^ in Vij (equation 
(2.18)) where boundary conditions (2.19) are replaced, respectively, by: 
1) the shear motion of the fluid between two neighboring inclusions: 



1 1 

Ui\dB* =9i = 2 e i> u AdBi = g 3 = ~2 ei > ( 2 - 21 ) 



2) the squeeze motion of the fluid between neighbors: 

u 2 \dB* = 9i = ^e 2 , u 2 \d B i = 9j = -\ e ^ ( 2 - 22 ) 

3) the rotation in the same directions: 

u 3 \dB' =9t = \{n\e 2 - n 2 e x ), u s \ dB j = g 3 = hn\e 2 + n\e x ), (2.23) 

4) the rotation in the opposite directions: 
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Fig. 2.10. The Poiseuille's microflow between motionless inclusions 



1, , , J 

2 l 



u4dB> =9 t = 77 W e 2 - «-2 e i), U4,\ dBj = g rj = --(n\e 2 + n\e x ). (2.24) 



Let u be the minimizer of (2.4) then 

u p := u- ui ub (2.25) 

minimizes Wu i} in the class 

V p = {v e H 1 ^) : V • v = in ZTy, v = on dB\ v = on dB j , 



Jlj ( v + u iub) ■ rids = (3* 3 1 



(2.26) 

The vector field u p describes the flow between two motionless inclusions, 
hence, it is natural to call it the Poiseuille microflow. 

Denote by (3* = {/3*j}iGiuBjeM and introduce the following set of discrete 
variables 

TV = {(U,w,/3*) : U\J satisfying (2.14) for % e B 

and (U,/3*) satisfying (2.17)} . ' 

It is straightforward to show that (see Lemma 4.2) 

W n = min V TD3n.Wn.A-), (2.28) 

where Vij is defined by 
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Vij = {v e H 1 ^) : V • v = in n ih v =U l + Ruj 3 {n[e 2 -n\ei) on dB i 
v = U j +Rcj i (n{e 2 - n{e{) on dB 1 (j e I) or v = f on dB j (j e B) 



1 

R 



v ■ nds 



(2.29) 

with U l , LO % and /3*j to be components of some no longer arbitrary but fixed 
(U, u>,/3*) G "72.*. Functions from Vy are defined in a single neck 27^ e 77. 
By direct computations vector fields 

«* = ^ + tn 

u sh = [(IP - U j ) ■ p l3 ] Ul + R(J + u j )u 3 , (2.30) 
w S9 = [(IT - U j ) ■ q ij ]u 2 , 

U per = R(u) 1 — W 3 )Ui + U p , 

are orthogonal with respect to the scalar product induced by the dissipation 
functional : 

Wn i} (U) = Wl7„ (U ah ) + Wily (U ag ) + Wj7„ (Uper), Wily ("0 = 0, 

(2.31) 

U = U t + U sh + U sq + U per . 

Physically, the decomposition (2.31) corresponds to three well-known types 
of microflows. Namely, (i) the shear type arises when a pair of inclusions 
cither rotates in the same direction or disks move into opposite directions 
(Fig. 2.8a, 2.7a), (ii) the squeeze type, when two inclusions in a pair move 
towards or away each other in thin gaps (Fig. 2.76), (Hi) and permeation of the 
fluid through the thin gaps between neighbors due to Poiscuille flow between 
motionless inclusions or rotation of neighbors in the opposite directions (Fig. 
2.10, 2.86). 

We are now ready to introduce a quadratic form which determines the 
discrete dissipation rate. To this end from now on instead of (3*j as a discrete 
variable we will use: 



flu - Pa - % W l + ■ p«] - ry - 



4i? 



(2.32) 



The reason for this replacement is that /3y is invariant with respect to Galilean 
transformation whereas f3*j is not (see remark in the end of Subsection 1). For 
example, if a constant vector U° is added to both U l and U j then the total 
flux /3*j changes while stays the same. Also, (3*j is the total flux through 
£ij of the entire flow (including the parallel translation, shear, squeeze and 
permeation) whereas is the flux due to the Poiscuille microflow solely. 
Finally, the use of simplifies the discrete dissipation form (equation (2.36)). 

The use of faj instead of leads to the replacement of the class TZ* (2.27) 
by TZ defined as follows: 
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7£={(U,u>,/3) : U\ lj 1 satisfying boundary condition on d(2 (2.14) for i e B, 
(U, (3) satisfying weak incompresibility condition(2.17), (2.32)} . 

(2.33) 

Introduce the overall discrete dissipation rate: 

J:=Q(U,w,^) = (u min^Q(U,w,/3), (2.34) 

where 

0(U,«,/J) = j;^ e , (2.35) 

iei jeAfi 

QijiU^U^w^LO^pij) = [(IT - U j ) -p ij + Rw i + Rwi] 2 C t ?8- 1 ' 2 
+ [(IT - U j ) ■ q ij ] 2 (Cp-V 2 + C^- 1 / 2 ) 
+ /3 2 (C l !5- 5 ' 2 + C^'r 3 / 2 + C^S- 1 ' 2 ) (2-36) 
+ R(u* - (Cp- 3 / 2 + C^- 1 / 2 ) 

+ R 2 {uj 1 - LjifC^S- 1 ' 2 , for j eJVifll, 

Q^IT, w \ /) = /3 2 [^<T 5/2 + Bp- 3 / 2 + Bp- 1 ' 2 

+ [(u l - u J ) ■ P ij + R^] 2 Bp- 1/2 

+ R 2 (u l -LjifBp- 1 / 2 

+ [(U* - U r ) ■ q ij ] 2 (Bp- 3 / 2 + Bp- 112 ) 

+ fa [(IT - U j ) ■ p« + iiu/] {Bp- 3 ' 2 + B^- 1 / 2 ) 

+ - ^){B%8- 3 ' 2 + S^tT 1 / 2 ) 

+ [(IT - CP) • + ifc;*] i?(u/ - c^')Biy~ 1/2 

+ [(£T - U j ) ■ q ij ] RaB l ^5- 1/2 , for je^.n B, 

(2.37) 

called the discrete dissipation rates, with coefficients C l £ , k = 1, ...,9, B^, 

R 

m= 1, . . . , 14, which depend on fj,, the ratio — and explicitly given by (A.l) 

d^ 

in Appendix 1. 

The solution of the discrete problem (2.34) is a set of discrete variables 
(U, u>,/3) e 1Z, where U represents the set of translational velocities of in- 
clusions, Q the set of angular velocities and (3 characterizes the Poiseuille 
microti ow in necks between neighboring inclusions. 
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Remark 2.2. Q(U, u>, (3), defined by (2.35), is a positive definite quadratic form 
(see Appendix 1). 

Remark 2.3. The agreement of the coefficients in (2.36) (explicitly given in 
(A.l)) with the previous results of [8] is as follows. Only coefficients C 1 /, C l 2 
coincide with the corresponding coefficients C*^, C l J p in [8]. This is because 
the coefficients in [8] are obtained by using the approximation of circular 
surfaces of inclusions by parabolas whereas in this paper we use the true 
circular surfaces. The main objective of [8] was capturing the strong blow up 
term of order <5~ 3 / 2 only, which requires coefficients C % j (C l $ ). The parabolic 
approximation does not bring any discrepancy in C l sp whereas it may bring 
a discrepancy in some other coefficients. Also since in [8] only the leading 
term was considered under special boundary conditions there was no need to 
consider the Poiseuille microflow. 

Both Theorems 2.1 and 2.2 are based on the following technical proposi- 
tion. 

Proposition 2.1. Suppose flp satisfies the close packing condition. Let the 
triple (U, oj, /3) € TZ solve the discrete problem (2.34)-(2.37). Then the follow- 
ing estimate holds: 

\W-1\ < M E C^+CzlU -U 3 \ 2 + C 3 ]T • (2-38) 

where Ck, k = 1,2,3, are dimensionless constants. 

The quadratic form Q defined by (2.35)-(2.37) can be written in the fol- 
lowing form: 

Q = Q i Jh + Qt q + Q i per + Q b p er+Q%, (2-39) 

where 

« = E E [(U'-U^-pV +Roj 1 + R^] 2 ^5- 1 '\ 

oi« = E E [(u i -ui)- q ^ 2 (cp-^+cp-^), 

+ R 2 (cv*-u,i) 2 Cp- 1 / 2 } 
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Q b P 



per 



= E E {fij [»i^ _5/2 + Bp- 3 / 2 + Bp- 1 / 2 



+Bl i p ij [(IT - U j ) ■ p ij + Ru l ] 5~ 3 ' 2 + B%l3 t3 R{^ - ^)<T 3 / 2 



+B§''/3y [(IT - U r ) ■ p 13 + Ru l ] S- 1 / 2 + B\\p l3 R(u l - lo^S- 1 ' 2 



+B%p iJ Rw t 5- 1 ' 2 + B l i [(IP -U j )- p 11 + Rcj'] 2 S' 1 ' 2 



+B% [(U l - U r ) ■ p 13 + Ruj 1 ] R(uj 1 - uj3)6- 1/2 



+B i iR 2 {J-^) 2 5- 1 ' 2 }, 



sq 



= E E {[(U i -U i )-q ii ] 2 (B$8- 3 ' 2 + tit's- 1 ' 2 ) 



(2.40) 



iei jeA/^nB 



+B% [{U l -U J )- q ij ] RaS- 1 ' 2 } . 



Remark 2.4-. Both Q™ r and Q b per are quadratic forms (see Appendix I). 

Remark 2.5. The discrete dissipation rate I makes physics transparent. We 
can see how microflows enter this dissipation rate. Indeed, the discrete dissi- 
pation form Q is presented as a sum of three motions corresponding to the 
decomposition (2.31). The decomposition (2.39) reflects the physics of the 
problem. Namely, for the interior necks the first term Q % s \ corresponds to the 
dissipation rate Wn^ (u s h) in (2.31) and describes the shear flow between in- 
clusions. This type of flow is produced by two motions: rotation in the same 
direction (Fig. 2.8a) and relative shear (Fig. 2.7a). The second term cor- 
responds to the dissipation rate Wn i:j (u sq ) in (2.31) and describes the local 
flow due to the squeeze motion (Fig. 2.7b). The third term Q™ er corresponds to 
the dissipation rate Wn id {u per ) in (2.31) and describes the "permeation" type 
flow due to the Poiseuille microflow between motionless inclusions (Fig. 2.10) 
and the rotation into the opposite directions (Fig. 2.86). Finally, the last term 
Qb describes the local flow in the thin gap between a disk and the external 
boundary. Here the decomposition into the parts corresponding to permeation 
Qp er and squeeze Q b in boundary necks is similar to one in interior necks. 

Remark 2.6. The asymptotics (2.38) is given in terms of the discrete dissipa- 
tion rates Qij . These dissipation rates are quadratic forms of the key physical 
parameters U, uj,/3. In particular, the overall discrete dissipation rate reveals 
the functional dependance of the overall properties of suspensions (dissipation 
rate, effective viscosity, effective permeability) on the microflows (the shear, 
squeeze and permeation between neighbors). 

3 Anomalous Rate of Blow up of the Dissipation Rate, 
Qualitative Conclusions, Discussions and Open Problems 



In this section we show that in 3D the Poiseuille microflow docs not occur 
(the fluid simply flows around the neck). In contrast, in 2D incompressible 
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fluid may need to permeate through necks. Necks separate fluid domain into 
disconnected regions, triangles, which may have different pressures pi, p 2 (Fig. 
3.1). In this sense, the 2D problem becomes more complicated than an anal- 
ogous 3D problem. Moreover, this type of flow in experimental 2D settings 
may lead to an observable physical effect. 




Fig. 3.1. Local flux due to pressure drop 



The techniques used in the previous study [8] of 2D flow were capable of 
capturing only the strong blow up term of order 0(S~ 3 ^ 2 ) of the asymptotics 
of the effective viscosity. The subsequent study of [10] reveals the significance 
of the weak blow up of order 0{6^ x / 2 ) in 2D. It was shown that it becomes the 
leading term in the asymptotics of the shear effective viscosity. The objectives 
of [10] was to study the degeneracy of the strong blow up term and evalu- 
ate the order of the magnitude of the next term which was shown to exhibit 
the weak blow up. The qualitative conclusion of the analysis of [10] was that 
the shear viscosity exhibits the weak blow up in both 2D and 3D. While this 
study highlighted the significance of the weak blow up term, the techniques 
of [8] did not allow to calculate it. As shown in Theorem 2.1 above the ficti- 
tious fluid approach suggested in this paper allows to calculate any singular 
term, and, in particular, the weak blow up. While a generalization of tech- 
niques of [8] for 3D suspensions leads to increasing technical difficulties, the 
fictitious fluid approach provides an appropriate tool to attack this problem. 
We anticipate that this approach would also be useful in a variety of similar 
physical problems (e.g. rigid inclusions in an elastic medium in both 2D and 
3D). Moreover, the analysis by this approach reveals the significance of the 
Poiseuille microflows. In this section we present an example which illustrates 
the following. For the suspensions of free inclusions the Poiseuille microflows 
contribute to the weak blow up. However, if an external field, which "clamps" 
inclusions, is imposed on inclusions then the Poiseuille microflow may result 
in a new type singular behavior of viscous dissipation rate (superstrong blow 
up). We also present an example of one disk in the fluid that can be clamped 
by the fluid flow with no external field. However, it is not clear whether this 
example can be generalized to an ensemble of inclusions. Our example may 
suggest that on a suspension of free inclusions the superstrong blow up may 
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occur only due to boundary layer effects and therefore for the large number 
of inclusions becomes negligible. 

1 Suspensions in a Pinning Field 

Example 3.1. The techniques described above can be applied to the problems 
of suspensions of non-neutrally buoyant rigid inclusions defined by 

(a) pAu = Vp, x e flp 

(b) V • u = 0, x e Of 

(c) u = U l + Ruj l {n[e 2 - n l 2 ei), x e dB\ i=l...N 

(3.1) 



(d) / o- 

Job* 



(u)n l ds = mg, i = l...N 

(e) I ri x er{u)n l ds — 0, i = l...JV 
(/)« = /, a; G 9f2 

where the domain is depicted in Fig. 3.2. This problem corresponds to 

JV 

minimization of Wr2 F («) + ^ mg • . The additional term docs not change 

i=l 

the analysis. 

Here we suppose that the density of the solid inclusions is p s and the 
fluid density is pf. Then the force on the inclusion in the left hand side of 
(3.1c?) counteracts the external gravitational field —mg, where g = (0,g) is 
acceleration due to the gravity and m = nR (p s — Pf) is the excess mass. 
We choose the external boundary condition / and the force exerted by the 
heavy disks on the fluid so that the inclusions do not move and fluid is forced 
to permeate through the thin gaps between motionless inclusions. The force 
exerted by the disks is equal to their weight TrR 2 (p s — pf)g, where we choose 
inclusions to be superheavy: 

p s = C(T 5/2 , for some C = C(p, R) > 0, (3.2) 

and the applied boundary data / is chosen so that 

'= (*£,)) GHl,2{dn) - (3 - 3 > 

where is some periodic function of (x, y) € dfl. Such a boundary data / and 
the asymptotics of p s are selected so that inclusions do not move (the gravity 
balances the viscous force). This balance can be found by solving an auxiliary 
problem similar to the one considered in [32] (for details see Appendix 4). 
Then the following proposition holds. 
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Fig. 3.2. Example of a domain occupied by a suspension in a presence of a pinning 
field (gravity) 

Proposition 3.1 (Superstrong Blow Up due to a Pinning Field). Let 

W be the overall viscous dissipation rate of the problem (3.1). There exist 
p s and f of the form (3.2) and (3.3), respectively, such that the following 
asymptotic representation holds: 

W = N (dS-V 2 +C 2 5- 3/2 +C 3 5- 1/2 ^) +0(1), as 5^0, 

where d = %i? 5 / 2 , C 2 = — n^R 3 ' 2 , C 2 = ^^-irpR 1 / 2 , 
4 ^ ' 160 P ' 17920 ^ 

and N is the total number of inclusions. 

The proof of this proposition is also given in Appendix 4). 

2 No Singularity of the Dissipation Rate due to the Poiseuille 
Microflow in 3D 

Example 3.2. In order to explain why the the Poiseuille microflow does not 
contribute to the singular behavior of the dissipation rate in 3D let us consider 
what happens with the fluid between two neighboring inclusions. Let K = 
[-L, L] 3 be a cube. The parts of two neighboring inclusions are modeled by 
the hemispheres attached to the top and bottom faces of the cube as shown in 
Fig. 3.3. Consider a 3D analog of the problem (2.1) with boundary condition / 
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to be a given constant vector V on two opposite sides of the lateral boundary 
and zero vector on the rest of the boundary (Fig. 3.3). Since the effective 
viscous dissipation rate W is bounded from above by the dissipation rate 
W(w) for any test function w it suffices to find w such that W(w) — 0(1). 



U=0 r 



u=V 




u=V 



Fig. 3.3. Poiseuille microflow in 3D 



Consider a "hourglass domain" W inside of the box containing two hemi- 
spheres as in Fig. 3.3. We choose the trial function w to be: 



0, when x E & 
W, when x G K \ 9 



(3.4) 



where W solves the Stokes problem with Dirichlet boundary conditions: W = 
on the boundary of ^ and W = f on dK (such a solution exists and is 
unique) . 

Evaluating the dissipation rate on such a trial function we obtain: 



W{w)<\\wf m{K) <C\\Vf H 



0(1). 



3 Boundary Layer Effects Leading to Superstrong Blow Up 

Example 3.3. Consider the domain fl = (—1, l) 2 which contains only one in- 

4 

elusion B. Decompose the domain outside B into necks 77^ and squares 



i=l 



[J Di as in Fig. 3.4(a). Choose 



! = 1 
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-1 N 



1 



-1 



0. 



Ci, ao 1 ndn=:r 1 



C 2 , da 2 n an =■. r 2 



C 3 , 9n 3 n9^=:r 3 



C 4 , 9Q4 n 9/2 =: A 



elsewhere 



where i = 1, . . . ,4 are smooth function having a compact support outside 
of such that Ci = 1 in A- Moreover, they satisfy the following symmetry 
condition: 

&(x,y) = y) = Caf-z, -y) = &(x,-y) and Ci(x,y) = d(y,x). 




1 X 



Fig. 3.4. One inclusion example: Boundary layer lead to superstrong blow up of 
viscous dissipation rate 



The fluxes through the parts of the boundary J\ and r 3 which are equal 
and we denote them as 

A) = |(-l + d), d = S + R(l-^=). 



Then the fluxes through I2 and Z4 are — /?o- 

Due to symmetry of the problem the inclusion does not rotate, that is, 
lo = and fa = fa = -fa = -fa (see Fig. 3.4(6)). Then 
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1 = min]T AS-^ 2 [U ■ qf + B^ X I 2 \U ■ pf + Ci? 2 <r 5 / 2 /3 2 
-Pi(fa+0i -04 + U!- U 2 ) - P2 {-/3 + p 2 - Pi - E/i - U 2 ) 
-ps(A> + A - ft - C/i + t/ 2 ) - p 4 (-A) + 0* ~ fa + E/i + E/ 2 ) 
- mjn{2AT 3 / 2 [C7? + L7f] + 2B ( T 1 / 2 [[/ 1 2 + + 40R 2 <T 5 / 2 /? 2 
-Pi(2A,+4/3 1 ) -^(-2/30-4/3!)}, 

where Pi, i = 1, . . . , 4 are the Lagrange multipliers corresponding to the weak 
incompressibility condition (2.17) and p\ = P3, p 2 — pa- 

Solving the Euler-Lagrange equations for this minimization problem we 
obtain 

u 1 = u 2 = o, Pi = -Y = li il - d) > 

which provides the following asymptotics: 

J = C(l — d) 2 ^ 5/2 . 

Thus, we see that the superstrong blow up can occur due to the boundary 
layer where Poiseuille flow is significant. 

4 Free Suspensions (no external field) 

In the above example 3.1 we demonstrated that the Poiseuille microflow dom- 
inates the asymptotics of the overall viscous dissipation rate and may cause 
the superstrong blow up of order 0(5~ 5//2 ). The key ingredient of this example 
is the presence of a strong external pinning field which "clamps" inclusions or 
alternatively the presence of the boundary layer as in Example 3.3. In typi- 
cal suspensions, with no external field, the inclusions are free to move. Then 
we expect that f3 is asymptotically small, as S — > 0, so that the Poiseuille 
microflow does not contribute to the superstrong blow up by Theorem 2.1. 
This observation is also supported by the analysis of the periodicity cell prob- 
lem with five inclusions in [9] . Below we present an example of a suspension 
with a hexagonal periodic array of inclusions and prove that for the exten- 
sional external boundary conditions the viscous dissipation rate exhibits the 
strong blow up (3(<5~ 3 / 2 ) since all /3y = 0. Such a hexagonal array in 2D is a 
representative of a densely packed array of disks. 

Example 3.4- Consider a square domain f2 = (— M, M) 2 , M > with hexago- 
nal array of disks B l centered at (x i} yi) G Q. We consider a periodic hexagonal 
array of such inclusions, as depicted on Figure 3.5). 

The boundary value problem (2.1) is supplemented with boundary condi- 
tions 

/ = ( _f ) on dn. (3.5) 
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dfl 




\ \ f t 7 

Fig. 3.5. Periodic domain occupied by a suspension under the extensional boundary 
conditions that exhibits the strong blow up 



We first prove that the overall dissipation rate Whex = 0(6 3 / 2 ) because we 
construct a trial vector field v such that Wq f (v) = 0(5~ 3 / 2 ). We thus obtain 
an upper bound for Whex ■ 

The construction is as follows. For each inclusion B z , centered at Xi = 
{xi,Vi) we prescribe the translational velocity to be exactly 



U' 



-Hi 



(3.6) 



and the rotational velocity uj 1 = 0. For such velocities the zero flux constraint 
(2.17) for every fluid region Aijk (see the right part of Fig. 3.5) takes the form 



Pij + Pjk + = 0, 
because elementary computations show that for every i, j, k 
(U l + U ] )p l] + (U j + U k )p> k + {U k + U l )p kl 



(3.7) 



0. 



Hence, we can choose all ftij to be identically zero. Then Theorem 2.1 implies 
that there exists a trial vector field v such that 

W Qp (v) =0(6-V 2 ), 

and, therefore, Whex = 0(<5~ 3 / 2 ). 

To show that for this array of inclusions Whex > C8~ 3 ^ 2 we consider a 
chain of disks that connects upper and lower boundaries of the domain tt: 
dfi + = {{x,y) : y — M} and df2~{(x,y) : y — — M}, respectively (see Fig. 
3.6). We choose this chain so that the y-axis of the coordinate system passes 
through the centers of disks in this chain. Then 

W = minQ(U,u;,/3) > (C 2 (T 3/2 + C 3 <T 1/2 ) niin ^ (U l 2 - U{f . 



chain 



Anomalous blow-up of the Viscous Dissipation Rate 



31 



I * * 

ogo 

°o° 

°o° 
ogo 

\ t t 



an." 



O 

0= 

o 
o 



Fig. 3.6. Chain of disks connecting the upper and lower boundaries of Q 



where the last minimum is taken over the disks in the chain (shadowed disks 
in Fig. 3.6). Then 

lin (U2 ~ T4f = A > 0, 



mm 

u 



therefore Whex > CS 3 / 2 . This leads to a conclusion that the overall viscous 
dissipation rate W per exhibits the strong blow up of order 5~ 3/>2 . 

It is known that for local movements of pairs of inclusions the squeeze 
type motion (Fig. 2.7b) provides the strongest singularity (of order <5 -3 / 2 ) 
whereas all other types of motions: rotations and shear provide a weaker 
singularity (of order <5 -1 ' 2 ). Thus, it is natural to expect that in a suspension 
where inclusions are free to move with the extensional boundary conditions the 
viscous dissipation exhibits the singularity of order <5~ 3 / 2 . The above example 
shows that this is indeed the case, that is, the superstrong blow up does not 
occur and the anomalous rate of 0{6~ 5 / 2 ) can be achieved if an external field 
is applied. 

There is also one more case where the superstrong blow up can be ob- 
tained. Namely, Example 3.3 shows that it occurs in the boundary layers. 
Hence the above 0{8~ z / 2 ) conclusion applies to bulk effective properties of 
free suspensions of large number of inclusions, when the boundary effects are 
negligible. 



4 The Fictitious Fluid Problem 

The use of the fictitious fluid approach immediately gives a lower estimate on 
the viscous dissipation rate W as follows. 

Lemma 4.1. ForWn defined by (2.10) andW defined by (2.3) the following 
inequality holds: 

W n < W. (4.1) 
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Proof. The idea of the proof is to observe that for any incompressible vector 
field u: Wn(u) < Wq f {u). Therefore, consider the minimizer u of (2.4) and 

W = /i [ D{u) : D{u) dx> f i I D{u) : D{u) dx. (4.2) 

since we disregard a part of the domain fi F . Hence, Wn(u) < W. 

We only need to show now that the restriction of the minimizer u on the set 
II is an admissible trial field for minimization problem of the fictitious fluid 
(2.10), that is, u\n € Vn- Indeed, u\n belongs to i? 1 (iT), is divergence-free 
in II and satisfies 

= / V • udx = I u ■ nds. 



Thus, 



min W n (v) < W n {u) < W = Wq f {u). 



□ 

Another significant advantage of the fictitious fluid problem is that the 
global minimization problem (2.10) can be split into two consecutive problems: 
one of them is on a single neck iZy , and the other one is a minimization 
problem on discrete variables (U, u),(3*) G TV . 

Lemma 4.2. (Iterative minimization lemma). Suppose Wn is defined by 
(2.10). Then 

W n = min V mmW nij (-). (4.3) 

Moreover, the minimizer of Wu i6 over Vij satisfies the following Euler- 
Lagrange equations: 

(a) ii Ait = Vp, x e llij , 

(b) V • u = 0, x e n^, 
(c') u = U i + Rw i (n\e 2 - n^), x e dB\ 
(c") u = U 3 ' + RLo j (n{e 2 - n{ei), x e dB^ , 

(d) ^-J u-nds = f3*j, 

(e) a(u)n = -p^n, xedll±, 

{ {f)u = f, xGdn^ndQ, 



(4.4) 



where p^ are the Lagrange multipliers for the weak incompressibility condition 



(2.17). 

Proof. Minimizing Wn(u) over Vn leads to the Euler-Lagrange equations 
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» 


/*Am = Vp, 


xen 


(b) 


V u = 0, 


x e n 


(c) 




-ni ei ), x e dB l , i=l 


(d) 


/ cr{u)n l ds = 


i = l...N 


JOB 1 




(e) 


/ n' x cr(u)n 1 ds 


= 0, i = l...N 


JdB* 




(/) 


/ m • nds = 0, 


i e I, j, fceJVj 


(.9) 


<r(u)n = p ljk n, 


a; e dAijk 


» 


u = f, 


x e 9/?, 



(4.5) 



where the pressure p(cc) arises from the incompressibility condition V • u = 

0, a; € II and "pressure constants" ptjk arise from weak incompressibility 
condition 

/ v ■ nds = 0. 

Given the boundary data / in (4.5) we uniquely determine (see Appendix 3) 
unknowns 

u, p, U\ uj\p llk , i e l,j,k e M- (4.6) 

Fix a neck 77^- and consider the problem (4.4) on it. For this pair of indices 

1, j take 



f,U\w\U>,uS,0* 



1 

R 



u ■ nds 



(4.7) 



found from (4.6). Using (4.7) as a given data, boundary value problem (4.4) 
can be solved uniquely (see Appendix 4). Due to the unique solvability of both 
(4.5) and (4.4) the pair (u,p) in (4.6) must solve (4.4) and 4 



Pi 3 =Pijk,P~i 



Pi jr. 



where triangles and Ay ra are adjacent to the neck 77^ (see Fig. 4.1). 

Hence (4.5) reduces to (4.4). 

At this point the completion of the proof would be trivial if we did not 
have /3*j. Indeed, in the entire domain ftp a result analogous to (4.3) is simply 



min Wq f (u) — min min Wq f {u) . 

(u,u>,u) ' (u,u>) y Ui when (u,w) fixed / 



Hence it only remains to show that for any given (U, u>) we can find at least one 
set (3* satisfying the weak incompressibility condition (2.17). To this end we fix 

4 For notational convenience we identify p+- = = pjjK = Pijh, and p+ = pj t . 
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(U, lj) and let u be the solution of the Stokes equation fiAu = Vp and V- u = 
in the domain Qp with the Dirichlet data u\ gB i = U l + RLL>' l (n\e 2 — n 2 ei), 

u\an = /■ Set /?*■ = ~ / u ■ nds. Hence, we obtain (3* = {/?*,-} such that 

(V,u>,f3*) £ TZ* . This completes the proof of lemma 4.2. 

Remark J^.l. For a given (U,u>) permeation constants (3* may not be unique. 
In fact, the choice of has N degrees of freedom where N is the number 
of inclusions. Indeed, (3* is found from solving a linear system A(3* = b 
where the number of unknowns equals the number P of interior necks and the 
number of equations equals the number of triangles, but there are only P — N 
linearly independent ones. Hence, the number of free parameters is equal to 
the number of inclusions. 





Fig. 4.1. (a) Pressures on the boundary of the triangle Aijt; (b) Two triangles 
adjacent to the neck Ilij 



Remark 4-2. The unknowns of the problem (4.5) are the velocity field u(x), 
the pressure p(x), the constant translational XJ % and angular io % velocities of 
the disk B l , i = 1, . . . , N and constants Pijk- Formally these constants appear 
as the Lagrange multipliers for the constraints (2.7). Similar to how the pres- 
sure p(x) appears as the Lagrange multiplier in the variational formulation 
corresponding to the Stokes equation, the weak incompressibility condition 
for the fictitious fluid, inherited from the original fluid, gives rise to the con- 
stant Lagrange multipliers, that one can regard as a constant pressure on 
the boundary of the fictitious fluid domain. This also motivates the notations 
Pijk- Thus, the fictitious fluid may be also interpreted as follows: an incom- 
pressible fluid occupies necks while triangular pockets are filled with "gas" of 
the constant pressure Pijk- Naturally, the unknowns of the problem (4.4) are 
the functions u(x), p(x) and constants pfp representing the velocity field, the 
pressure in the neck 77^ and the constant pressures on the lateral boundary 
dn^j, respectively. 
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Remark 4-3. A major difficulty in applying previous 'one-step' discretization 
techniques [8, 12, 13] to vectorial problems is the presence of integral con- 
straints in the dual variational formulation. In our 'two-step' discretization 
approach, due to the Iterative Minimization Lemma, the inner minimization 
problem has Dirichlet boundary conditions on inclusions and therefore neither 
this problem, nor its dual have any integral conditions. On the other hand, 
due to (4.3) second minimization implies that these integral conditions are 
automatically satisfied. 

Lemma 4.3. Suppose Qf satisfies the close packing condition. Then for Wn 
defined by (2.10) and W defined by (2.3) the following inequality holds: 



' (4 - 8) 

where (U, a),/3*) minimizes (4.3). 

The proof of this lemma relies on the technical construction that appears in 
Subsection 1 and therefore is put in Subsection 2. 

Remark 4-4- Lemma 4.3 is the only place where the close packing condition 
is necessary to obtain the desired estimate (4.8) because we needed uniform 
Lipschitz regularity of triangles A^. 

As a corollary of Lemmas 4.1 and 4.3 we have the main result of this 
section: the accuracy of approximation of the overall viscous dissipation rate 
by the dissipation rate of the fictitious fluid given in the following proposition. 

Proposition 4.1. Suppose Of satisfies the close packing condition. Then 

\W-W n \<»(j2 E C ^ R % + C *\U U j \ 2 + C 3 E R2 (^) 2 

(4 - 9) 

where (U, <D,/3*) minimizes (4.3). 



5 Discrete Network 

In the previous chapter we described a discrete network that arises from the 
fictitious fluid approach. Indeed, the equation (4.3) in view of (4.4) is 

W n = min V Wn.Au) =: min WOU,u,0*), (5.1) 

where u is the solution of (4.4) and W is a positive definite quadratic form 
of (U,«,/3*). 

Our next objective is to find coefficients of W asymptotically as charac- 
teristic distance S — > 0. We have the following result. 
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Lemma 5.1. Suppose Wn is defined by (2.10), X is defined by (2.34), W is 
defined by (5.1), and Q is defined by (2.33)-(2.37). Then 

|w(u,u,,/T)-Q(u,u,,/3)| < M ]T Ci&ul + wir - u*\ 2 + c 3 ]T r2 ( 

j&Afi ieluB 

(5.2) 

for any (U, u>,/3*) € 7?.* and /3y related to (3*j through (2.32). In particular, 
\W n -l\<JY J Y. C i R %+ C >\U i -U j \ 2 + Cz E R2 &? 

' (5 - 3) 

where (U, o),/3) is the minimizer of Q. 

This Lemma shows that coefficients of W(U, u>, /3*) tend to infinity as 
(5^0 because the corresponding coefficients of Q(V, w, (3) are asymptotically 
large and given in (2.36)-(2.37). The proof of this Lemma is given in Subsection 
1. 

Combining Proposition 4.1 and Lemma 5.1 we obtain the claim of Propo- 
sition 2.1. 

In order to prove Theorems 2.1 and 2.2 it remains to show that the error 
term of the right hand side of (5.3) becomes relatively small compared to the 
overall discrete dissipation rate I. In order to show that we prove in the next 
lemma that I — > oo as 5 — > 0. More specifically we have the following result. 

Lemma 5.2. Suppose Qf satisfies the close packing condition. Then there 
exists a constant C > such that for every {U 1 }: 

I>^^ dS-V^U* -U j \ 2 + C 2 5- 1 ' 2 R 2 {J + u>i) 2 

iei jeMi (5.4) 
+C 3 S- 1 / 2 R 2 (uj t -uji) 2 +C 4 S- 5/2 (3f v as 5^0. 

Proof. We write 

J = E Ctf-WftlT - U j ) ■ q l] ] 2 

+C 2 S- 1/2 [(U l - U j ) ■ p ij + Roj 1 + Rcji} 2 + CaS-WR 2 ^ - u 1 ) 2 
+C 4 5~ 5 / 2 R 2 f3 2 j . 

Using (a + b) 2 > 3/4a 2 - 3b 2 we have 

i>^^ Ci<r 3 / 2 [(cr - u r ) ■ q l] ] 2 - 3C 2 5- 1/2 \u l - u°\ 2 

+- A C 2 5- 1 ' 2 R 2 {J + <J) 2 + CsS-^R 2 ^ - ^) 2 + C 4 <r 5/2 i? 2 /3 2 . 
4 J 

By discrete Korn's inequality (Lemma A.l, Appendix 2) when 5^0: 



a/) 2 
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dS-^liU 1 - U j ) ■ q 11 ] 2 - 3C 2 8- 1 / 2 \U i - U 3 \ 2 > C5- 3 / 2 \U l - U 3 \ 2 . 
Hence, (5.4) holds. 

From Lemma 5.2 and Proposition 2.1 we have Theorem 2.2, and Theorem 
2.1 follows from Theorem 2.2 and Proposition 4.1. 

6 Proofs of Lemmas 5.1 and 4.3 
1 Proof of Lemma 5.1 

Boundary value problem on a neck. In the local coordinate system the 
neck Ilij (see Fig. 6.1) is given by 

Ilij = {x = (x,y) £ M 2 : < x < j±, -Hu(x)/2 <y< H l3 {x)/2) , 

(6-1) 

where \^f- — 7^-| is the neck width, and H i3 is the the distance between dBi 
and dBj defined by 

Hij(x) = Sij + 2R- 2y / R 2 ~x 2 . (6.2) 




Fig. 6.1. (a) The neck Ilij connecting two neighbors B 1 and B^ , (6) the distance 
H (x) between two neighbors, (c) the boundary dllij of the neck 



Fix (U, u>,/3*) £ TV. In the chosen coordinate system the boundary con- 
ditions on dB 1 and dB^ for functions defined in (2.30) are written as: 



38 Leonid Borlyand, Yuliya Gorb, and Alexei Novikov 



u.fclam = G\ - -{U{ - U{)e x + i?K + uP) [5V 1 , TP 



2R 



u sh \ dBi =G{ = -\(U{-U{)e 1 +R{u i +^)[ W 1 TP 



(6.3) 



2R 



u sq \ dB * = G\ = \{XJ\ - U()e 2 , 
u sq \ dB] =& 2 = - l -(U{-U{)e 2 , 



u per | dBi 



= G\ = R{uo l - uji) (V 1 - 3* j 

V 2R. ) 

V ~2R 



2R 

UperldBi = G J 3 = R(UJ 1 — LU j ) I 2 V ^ _ "R 7 



(6.4) 



(6.5) 



and u s h, u sq , u per are minimizers of Wn^ over the following sets, respectively: 
V ah = |u e ff 1 (77y) : V • u = in 77y, v = Gi on 55', r> = G{ on 

(6.6) 

F S9 = {»£ H 1 ^) :V -v = in n i:j , v = G\ on 95% u = G 2 on 5^} 

(6.7) 

Vp er = € H x {n i3 ) : V ■ v = in n i:j , ,v = G l 3 on dB\ v = Gf 3 on <95 J 



(6.8) 

Variational Duality. For any (U, oj,(3*) £ TZ* and any 77^ £ II we let 

Wy (IT*, ^tf) : = ™ in ^ (•)■ ( 6 ' 9 ) 

By orthogonality (2.31) 

WijilP, U j ,u\ ^,/3*-) = Wir„ (« ah ) + Wn tj (u sq ) + W nt] (u per ). 
Using variational duality (see e.g. [11, 12] ) 5 

W(u 8 h) = maxWi*(S), (6.10) 

W(u sq ) = maxW 2 *(5), (6.11) 

«SfE F 

W{u per ) = max W 3 *(S), (6.12) 

5 Since the construction is done in the neck ITij hereafter we drop such a subscript 
in the functional Wn i} = W. We also drop the indices i,j in Hij(x) = H(x). 
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where the dual functionals are defined by (k = 1, 2) 

(trS) 2 



W* k {S) = f Gl- Snds + f G{- Snds - ±- f 



S:S 
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dx, 



Wi 



'(S)= f G\- Snds + [ Gl- Snds 

JdB i JdBi 

R0T< 



4/U J n 



n ■ Snds 



S :S 



(trS) 2 



n ■ Snds 



(6.13) 



(6.14) 



dx 



over the sets 



F={se M. 2x2 ,S = S T , Sij G L 2 {n i3 ) : V • S = in 77^ Sri = on <977±} . 

(6.15) 

F per = {SG R 2x2 , 5 = S T , S tf G L 2 (i7y) : V • S = in 77^, 

Sn = X ±non 577±} , 

(6.16) 

and xtj an( i Xy are arbitrary constants. In (6.14) the numbers H(j^) are the 
lengths of the lateral boundaries, that is, 77(7^) = \dll^\, 77(7^) = |d77^|. 
It is straightforward to check that the maximizer of problem (6.10), (6.13), 
(6.15) is the stress tensor <r(u s h), the maximizer of (6.11), (6.13), (6.15) is 
er(u sq ) (see e.g. [8]) and the maximizer of (6.12), (6.14), (6.16) is cr(u per ) (see 
Appendix 5). 

Proof. By direct and dual variational principles for any S\ G F, S2 G F, 

S3 G F per and v\ G V s h, v 2 G V^ g , G V per we have the following bounds: 

WfOSi) + VK 2 *(<S 2 ) + W 3 * («S 3 ) < < + VK(t> 2 ) + VK(t> 3 ). (6.17) 

Therefore, if we find trial fields and Uj such that the following three in- 
equalities hold (with universal constants d, i — 1, 5) 

\W( Vl ) - W?(5i)| < n (Ci[([/ 1 - U j ) ■ + Ruj 1 + Ruji} 2 + C 2 R 2 (cu l + ^) 2 ) 

(6.18) 

\W(v 2 ) - W*(S 2 )\ < iiC 3 [(U l -U j )- q^] 2 =: E 2 (U\ U j ), (6.19) 

\W(v 3 ) - W;(S 3 )\ <^(C 4 R 2 (u; t -^) 2 +C 5 R 2 l3 2 3 ) =: E 3 (w\ w j , /%), 

(6.20) 

then for v = u t + v\ + v 2 + v 3 we have 

\W(u) - W(v)\ < \W( Vl ) - W7(5i)|+|W(t» 2 ) - W?(S 2 )\ + \W(v 3 ) - W?(S 3 )\ 
< ndtfj + C 2 \U l - W\ 2 + C 3 R 2 (uj 1 ) 2 + dR 2 ^) 2 , 
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and Lemma 5.1 will follow from summation of the above inequality over all 
necks and explicit evaluation of W(v\), W(v 2 ) and W(v 3 ). 

The construction of <S^ and Vi for the interior necks is done, for simplicity, 
for symmetric necks (7^ = 7^ =: 7»j) and it is as follows. 



t>i = (U{~U{+Ruj 1 + Ru) j ) 

yK(x) 



yG(x) 



(6.21) 



where functions F, G, K and M are found so that v\ satisfies (6.3) on 
dB\dBi, that is, 

(6.22) 



'1-fw-l 



#(*) ' 



G(s) - C 



5l= ^G(z)-G -yG'(x) J in ^ (6 ' 23) 

where the constant C is chosen so that S\ £ F, namely C = —j- — -. 

ti \"fij ) 

• 9 (*.y) = (^-^)(_^^t)) in ^' (6 - 24) 

where G(x) and F(x) are chosen so that v 2 to satisfy (6.4) on dB l , dB J , that 
is, 

3r 2r 
G ^ = -2^)' F ^ = ^)- ^ 

(S 2 )n = - C/|) ^G' + 9y 2 F' -qJ'f-C^- ^C 2 x 2 y 2 + C 3 + C 4 y'< 
(S 2 )i2 = (5 2 ) 2 i = M(^2 - ^2) (6y^ - V°G" - 3y 3 ^" + C x xy + C 2 xy 3 ) 

(S 2 ) 22 = U 3 2 ) (-3G' - 3y 2 F' - 6 |* F + ^G'" + ^ 4 F"' - G^ - G 2 

(6.26) 

where constants G, (« = 1..4) are chosen so that S 2 n = on 977^. 
„ P(x)+3y 2 0(x) \ f K(x)+3y 2 M(x) 

(6.27) 
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where functions P, Q, K and M are found so that v 3 satisfies (6.5) on 
8B\dBi, that is, 



P(x)-. 
K{x) 



3x 2 3H(x) 1 
ARH{x) + 16R + 2 



0te) = — — 

i? 2 ' WW iJff3(a.) 4flff(a;)' 



2iJ(x) ' 



M (jc) = - 



ff 3 (x)' 



(«S 3 )n = »(p' + 9y 2 Q' -gJq-C 

(«S 3 )i2 - (5 3 ) 2 i - M (6yQ - |/P" - 3y 3 Q" + C lX y + C 2 xy 3 ) , 



(6.28) 



(s 3 ) 22 = n ( -sp' 3y 2 g' - 6 / q + y -p"' + ^V" -c^-c 2 y -\, 



where 



4 

(6.29) 



G(x) = R{lo 1 - oj j )P(x) + /3ijRK(x), F(x) = R(uj 1 - lu 3 )Q(x) + j3 i:j RM(x), 

(6.30) 

and constants (i = 1..4) are chosen so that <S 3 n = x ±n on f° r some 

constants x^- 

In Appendix 1 we verify that the chosen trial fields indeed satisfy inequal- 
ities (6.18), (6.19), (6.20), and also show that 

W{vi) = ^^R 1/2 (U{ - U{ + Rj + RtJ)% 1/2 

+C^(Ul - U{ + Rlu 1 + Rw j ) 2 + C 2j ui? 2 (w 4 + u j ) 2 (6-31) 
+C 3 fi(Ul - U{ + Ruj i + Ruj J )R 2 {uj l + 



W(v 2 ) = (Ul - Uif 



and 



1/2 



320 



+ C^(Ui-Ui) 2 , 
(6.32) 



+i? 2 K-^") 2 ^ M i? 1 /2^V2 

+C 5 i? 2 K - J) 2 + C 6 i?V - uP)Pa + C 7 i? 2 /3 2 . 



(6.33) 

Coefficients of the quadratic forms Qij (2.36) can be read off from the 
formulas (6.31)- (6.33). Construction for the boundary necks and computation 
of the coefficients of (2.37) is identical and it is given in Appendix 2. 
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Remark 6.1. Note that the first term of the decomposition (6.27) indicates the 
contribution from the fluid motion due to the rotation of neighboring disks in 
opposite directions while the second term describes the velocity of the fluid 
between motionless inclusions that is the Poiseuille microflow in a neck. Con- 
stants )3ij , called permeation constants, characterizing this microflow, are the 
local fluxes through the line segment tij connecting two motionless inclusions 
and given by (2.32). 

2 Proof of Lemma 4.3 

In the previous subsection we showed that for any (U, u>,/3*) € TV there 
exists a trial vector field v n on the fictitious fluid domain such that 

Vo\dB* = U l + Ruj l (n\e 2 - n?,ei), v \ 9n = f, 

and for any neck 77^- 

< W nij (v )-Wn tj (u) < C 1 \U i -U j \ 2 +C 2 R 2 (u i ) 2 +C 3 R 2 (u j ) 2 +C 4 R 2 l3 2 j , 

(6.34) 

where u is the minimizer of (2.10). 

Hence, in order to prove Lemma 5.1 it suffices to show that we can extend 
t>o in every triangle Aij k continuously so that this extension, called w, satisfies 

W Al]k (w) < d\U* - u j \ 2 + c 2 \u j - u k \ 2 + c 3 \u k - c/f 

+C a R 2 (uj 1 ) 2 + C 5 R 2 (cui) 2 + C 6 R 2 (uj k ) 2 
+C 7 R 2 pl + C 8 R 2 p 2 k + C 9 R 2 (3l, 

in every triangle Ay fc . We construct this extension by assuming that w solves 
the Stokes equation on A^k 6 : 

{(a) fiAw = Vp, in A ijk 
(b) V • w = 0, in A ljk (6.35) 
(c) w = v , on dA ijk 

For (6.35) by [26] we have 

W^ k <C\\v -U\\ 2 m/2(dAijk) , (6.36) 

where U is any constant vector. The constant C may in principle depend 
on the Lipschitz constant of A^ k , but, since the close packing condition is 
assumed, angles in all Aij k are bounded uniformly from and ir, and therefore, 
C in (6.36) is universal. 

Let U = l/'S(U l + U 3 + U k ). Then (v - U) on dA ijk is explicitly given. 
Namely, 

6 Here we discuss interior triangles, estimates for the triangles at the boundary are 
identical. 
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(t» - U)\ Bn + = l -(U l - U k ) + l -(U ] - U k ) + [(IT - W) ■ p« + Rj + Rj] v\ 
+ R(uj 1 + ujJ)vj + [{U l - U j ) ■ q ij ] v 2 + R(w l - u j )vl + f3, 3 v 2 , 

where dA i3k = dll^ U dllf k U dll^ (see Fig. 6.2) and v\, vf, v 2 , v\, v 2 are 
polynomial vector fields independent of (U, w,/3*) and 5, and given in (6.21), 
(6.24), (6.27), and f3. t] is defined by (B.16). Therefore, 




Fig. 6.2. Boundary of the triangle Ayt 



\\v Q u\\l 1/2(dAijk) < c^u* U>\ 2 + C 2 \U* -U k \ 2 + C 3 \U k U'f 

+C i R 2 {u i f + C 5 i?V') 2 + C 6 R 2 {u> k ) 2 
+C7R 2 /3 2 j + CgR 2 f3 2 k + CgR f3 ki . 

(6.37) 

Combining (6.37) with (6.36) we obtain 

E E ^ iJt W<EE C!\U* - Ui\ 2 + C 2 \Ui - u k \ 2 + c 3 \u k - u*\ 2 

+C 4 R 2 (uj') 2 + C 5 R 2 (lu^ 2 + C 6 R 2 (tu k ) 2 
+C7.R /3y + C 8 R 2 (3 2 k + CgR 2 (3^. 

Finally, the above inequality and (6.34) give the desired result (4.8). 
7 Conclusions 

In this paper the asymptotic formula for the overall viscous dissipation rate 
is obtained, where all singular terms are derived and justified. This is done 
by developing a new technical tool - the two-step fictitious fluid approach. 
Such approach is expected to be helpful in evaluation of effective properties 
of various highly packed particulate composites. The obtained asymptotics 
provides for a complete picture of microflows, while previous studies, men- 
tioned in Introduction, gave only partial analysis of microflows and singular 
terms. In particular, a new term due to the Poiscuillc microflow, which was 
not taken into account previously, is obtained. It is shown that this Poiscuillc 
microflow does not contribute to singular behavior of viscous dissipation rate 
in 3D. In contrast, in 2D it may result in an anomalous rate of blow-up (of 
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order <5~ 5 / 2 ). Indeed, such a rate of blow-up is obtained in the presence of 
external field (e.g. gravity) or due to the boundary layer effects. Our analysis 
suggests that in absence of an external field the anomalous blow-up does not 
occur. 

The obtained asymptotics expresses the continuum dissipation rate in 
terms of a discrete dissipation rate, and the latter reveals dependence on 
the key physical parameters. 

We remark that this work leads to a somewhat surprising observation that 
suspensions are actually harder to analyze in 2D than in 3D. As we mentioned 
above, the Poiseuille type microflow is significant in 2D and it is negligible 
in 3D. The key reason here is topological: in 2D thin gaps between closely 
spaced inclusions partition the fluid domain into disconnected regions, which 
is not the case in 3D. Hence, in 2D permeation of fluid between two inclusions 
contributes into the singular behavior of the overall viscous dissipation rate. 

Finally, we note that 2D mathematical models were often used to analyze 
qualitative behavior of 3D problems in order to reduce the analytical and 
computational complexity of the problem. Our work clearly shows limits of 
validity of such modeling. 



The coefficients , k = 1, . . . , 9, and B]^, m = 1, . . . , 14 appearing in (2.36) 



A Appendix 

1 Coefficients of the quadratic form Q 

■>ij 

'k i 

and (2.37), respectively are given by 

c H^(£f' C H^GI) 3/2 ' c *=^{£ 1/2 

Mi 9 f R \ 5/2 Mi 99 / R \ 3/2 Mi 29241 / R \ 1/2 
ei = 4^{d- j ) ' ^leO^UJ ' 6 ? = 17920*" Uj 

^■=-^(|) 3/2 '^=^(^) 1/2 , <*=H£ X 1/2 
B - 4 ^(i) 1/2 ' B ?=Hi) 1/2 > *° =6 A£ 3/2 

/ R \ 3//2 3 / D V 1 / 2 / n \ V 2 

p \ 1/2 / p \ 1/2 



B\% = -3tt^ ( — J , S J / 4 = 67r/i ( -p 



(A.l) 
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2 Discrete Korn Inequality 

Lemma A.l (Discrete Korn's Inequality). Suppose Qp satisfies the close 
packing condition. Given {U^ \i € B} there exists a constant C > such that 
for every {U\i G 1} the following inequality holds: 

E E [(^ - u °) ■ > C EE^- u3 ) 2 - ( A - 2 ) 

iei jeTVi iei jeM 

Proof. By contradiction we assume that there exists {t^^} such that 

E E [(^ - u i) ■ ^ < I E E - u{f. 

iei jeAfi iei jeM 

Then we have two options here. 

Case 1 . There exists a constant M > such that for alH e E and k e N: 

Case 2 , There exists in e I and a subsequence {k n } such that C/^ — ► oo 

as fc„ — > oo. Without lost of generality we set: |{7l° | = max It/ 1 . |. 

n iei n 

First we prove the Case 1. By Bolzano- Weierstrass lemma there exists a 
subsequence {k n } and a set of vectors {V l ,i E 1} such that U\ n — > C/ 1 as 
fc„ — > oo for iei and for V 4 we have: 

(ai) (V* - F j ) • g« = 0, i.jelUB; 

Recall that in this paper we focus on the quasi-hexagonal array of inclu- 
sions which triangulates the domain Q and in which every disk typically has 
six neighbors. For such an array the condition (ai) implies that all disks, in- 
cluding quasi-disks, are moving as a rigid body, that is, translate and rotate 
as a whole. Hence, there exist numbers a, 0, <p such that 

^/aWcosy, sinyA /x\ (A 3) 

\PJ \-sm<pcos<pJ \yj 

However, (&i) implies that a = = 0, cos 95 = a = 0, b = — s = simp which 
are excluded by our boundary conditions (2.2). 

For the Case 2 denote by dk n := |C7^ | = max \U\ n | — » 00 as fc„ — > 00. 

Consider = Hk. We know that |v£j < 1, |V*?J = 1. Hence, there 

exists a subsequence {fc„} and a set of vectors { V 1 , i E 1} such that 
(a 2 ) Vfe n — > V* as fc„ — > 00; 

(62) = i; 

(c 2 ) V 4 = 0, for i e B; 

(<fc) (V - V*) ■ q ij = 0, i,j €IUB. 

As before, the condition (d 2 ) implies that all disks move as a rigid body, 
that is, (A.3) holds. Using (c 2 ) we have that a = = <p = 0, hence V 1 = 
for i e I U B which contradicts (62). 
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3 Construction of the network 

Here we describe how the network Q is constructed. 

We follow [12] where the notion of neighboring inclusions were defined 
via the Voronoi tessellation (diagram) [5] of the domain £2f (Fig.A.l). For 
an arbitrary distribution of the centers Xi of inclusions B l , called sites or 
vertices: 



X = {x t . e n, i = 1,...,N}, 
the Voronoi diagram is the partition of the domain Q into non-overlapping 



Voronoi cells V(xi). Each V(xi) is the set of all points in J? that are closer 
to Xi than to any other site from X. The situation when site Xk is located 
near one of boundaries of J? is treated similarly to the one in [12, 13]. Namely, 
one can consider the reflection x^ of x^ along this boundary and construct 
the corresponding Voronoi cell using the obtained auxiliary cite X)~" outside 
of Q (see Fig. A. 2, where the hatched region is a neck), introducing a new site 
xy on the boundary, which is the center of a line segment B k in dSl, called 
a quasidisk (see also [12, 13]). Two disks are called neighbors if their Voronoi 




Fig. A.l. Voronoi diagram and Voronoi cell of the site 




Fig. A. 2. The quasidisk B k centered at x k i on the boundary of fl 



cells share a common edge. 
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Fig. A. 3. A "triangle" near the boundary of Q 

Now we can define a network Q as introduced in Subsection 2. 

This graph is simply the Delaunay triangulation of X. Recall that the 
geometric dual of the Voronoi diagram in two dimensions is a Delaunay tri- 
angulation [5]. The basic element (or cell) of the Delaunay triangulation is a 
triangle that may degenerate into quadrilateral or n-gon, when four or more 
disk centers lie on the same circle. 

The neck-triangle partition is a decomposition of the domain flp into sim- 
ple geometric objects: necks and triangles, based on the central projection. 




(a) (b) 

Fig. A. 4. (a) Central projection projects the Voronoi edge AB onto arcs A'B' and 
A"B" , (6) the neck ily 

The central projection partitions the boundary of a disk using the cor- 
responding Voronoi cells as follows Hi : dV(xi) — > dB 1 , i € I, where 

Ki(x) = R(X - Xi)/\X - Xi\. 
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As shown on Fig. A. 4, the neck II fj connects arcs A'B' with the arc A" B" , 
obtained by the central projection of the same edge AB. We remark here that 
with this construction the edges of a neck are always parallel. Also a neck 
is necessarily symmetric with respect to the line connecting centers of the 
neighboring disks. The set of all necks in Qp is U (the shadowed region on 
Fig. 2.4), called the fictitious fluid domain. 

Once necks are defined the complementary part of the domain is the union 
of triangles (see e.g. Ayt on Fig. 2.4) and possibly trapezoids near the bound- 
ary dfl. As mentioned above, by slight abuse of terminology we call such 
trapezoids "triangles" . The set of all triang les is A = (J A ijk . 

Thus, we obtain the neck-triangle partition of the fluid domain: flp = 
A U II. Such a partition of the fluid domain Q F allows us to decompose the 
functional Wq f as follows: 

W = WnA') = w n{~) + WA(-) =: W n + W A . 




Fig. A. 5. (a) The neck Ilij connecting two neighbors E 1 and B 3 , (b) the distance 
H (x) between two neighbors, (c) the boundary dlJij of the neck 



4 Auxiliary Problem for (3.1) 

Here we prove Proposition 3.1. First we show that it is possible to choose such 
a density p s and boundary data / that disks do not move. 



Theorem A.l. There exist p s and f such that U % = and ui l — 0, i = 
l,...,iV. 
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Proof. On the periodicity cell Y = (—5, 5) x ("515) (Fig- A. 6) consider: 



' (a) pAv = Vp, 
0) V • v = 0, 
(c) v = 0, 

(e) v is y — periodic 



cc e y\s 

x e Y\B 
x edB 



udx 



(A.4) 



This problem is equivalent to the one considered in [32]. 

The existence of the unique solution to (A.4) can be verified by the classical 
energy methods. Observe that the solution v = (v\,V2) of (A.4) satisfies the 
equations (a), (b) of the problem (3.1), and the condition (3.1c) (with U l = 0, 
uj 1 = 0) for Q = Y (i.e. N = 1). Using symmetry of the problem one can show 




that v also satisfies the balance of the angular momentum (3.1e) with fl = Y: 

I cr(v)n x nds I = / ([cr(v)n]in2 ~ \cr{v)n\xn\) ds) =0, 
Job V JdB J 

for some a?M which is not, in general, zero (this integral is the force exerted 
by the fluid on the inclusion B). 

We estimate a in terms of 5. For this we compare problems (3.1) and (A.4) 
for u and v, respectively. If we choose / in (3.1e) to be a restriction of v on 
dfi and a = mg = TrR 2 (p s — pf)g then v solves (3.1). Hence, u — v in ft = Y, 
where, in particular, the translational and angular velocities are zeros: [7 = 0, 
u = 0. 

Consider now the problem (3.1) for u obtained by a periodic extension of 
v defined in Y to N cells. If we now choose p s and / as above, namely, 
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then uniqueness implies that all inclusions arc motionless (due to (A.4c)). 
□ 

Theorem A.l shows that if we choose the density p s and the external 
boundary conditions / in (3.1) so that inclusions are motionless then the 
only microflow between neighboring inclusions is permeation (or the Poiseuille 
microflow). We now show that the overall viscous dissipation rate of (3.1) 
exhibits a superstrong blow up. Hereinafter we use the equivalence of (3.1) 
and (A. 4) extended by periodicity to fi. Before we obtain the asymptotics of 
the dissipation rate W for this problem we find an asymptotics of p s as S — > 
and the form of /. 

Lemma A. 2. The following asymptotics holds for the force exerted by the 
fluid on the disks B defined by the right hand side of (A.J^d): 

a = C5~ 5/2 + 0(<T 3/2 ) with some constant C = C(m, R, N) > as 5 -> 0. 

Proof. Consider an auxiliary function w = v — (v) that solves: 

' (a) pAw = Vp, xe Y\B 

{b)V-w = 0, xe Y\B 

(c) w = , xedB (A.6) 

(d) (w) = 

(e) w is periodic 

Then if we multiply (A. 4a) by w and integrate by parts over Y\B we obtain: 
= — / D(w) : cr(v)dx + I w ■ (r(v)nds + I w ■ cr(v)nds. 

JY\B JdB JOY 

Taking into account that D(w) = D(v) and (A. 6c) we continue 

= — / D(v) : cr(v)dx + / [<T(v)n] 2 c?s + / w ■ cr(v)nds 

JY\B JdB JdY 

= — W y \b( v ) + a + / w ■ cr(v)nds. 

JdY 

Due to periodicity of the function w, indicated above properties of the com- 
ponents of the vector a(v)n and the symmetry of the domain Y we obtain 

that / w ■ cr(v)nds — 0. Hence, a — Wy\b( v )- Therefore, in order to ob- 

JdY 

tain an estimate for a with respect to S we have to study a problem of the 
asymptotics of the overall viscous dissipation rate of the suspension given by 
(3.1). Since inclusions in problem (3.1) are motionless then the fluid permeates 
through thin gaps between them. So in view of formula (2.39) of Theorem 2.1 
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the viscous dissipation rate of this suspension exhibits the superstrong blow 
up of order 0(<5~ 5 / 2 ) if we prove that /3y = O(l) for at least one pair 
i G l,j G 

Recall that / now is given by (A. 5) and v-y{— \,y) = vi(^,y) due to 
periodicity of v. But v% is odd in x, hence, v\ does not depend on x. Similarly 
v\(x,— |) = V\(x, 5) = due of oddness of V\ in y. Thus, V\{x, y) = 0. 
Moreover, due to the periodicity of V2(x, y) we have that the second component 

f2{x,y) of / is even in both x and y. Hence, / = (j, ^ where fa is an 

even function of both x and y on dil. 

Since for p s and / given by (A. 5) we have U % = 0, and uj 1 = for all i, the 
problem (3.1) becomes the problem of minimization of Wn F (v) introduced in 
(2.3) over the class V defined by (2.4), where U % = 0, and lo 1 = for all i, 
considered above. As mentioned above a = W asymptotics of which is given 
by Propostion 2.38. To find an asymptotics of W we need to know the flux 
Pij through the line £ij joining two neighboring disks B l and B 3 . In order to 
find this parameter consider a periodicity cell Yf as in Fig. A. 7(a). We note 
that due to divergence-free property and periodicity of it one can obtain that 
Pij is equal to the flux through any horizontal line £h connecting the lateral 
boundaries of Yp' 



/%4 



u 



nds = — 
R 



5/2 x 

u%{x, 0)dx = — 

-5/2 



R 



u 2 (x,y)dx, 



for any fixed y G (—1/2, 1/2). Using this fact, u = on dB, (u) = 



and 




B n 




^4 




I, 








Bt 








B m 


\ ( 





(a) 



Fig. A. 7. Periodicity cell Yf 



Fubini's theorem we have 
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1 = 7777/ u 2 (x,y)dx = / dy u 2 (x,y)dx = / dy u 2 (x,y)dx = Rfaj. 



In Fig. A. 7a £ij is shown as a horizontal line. There are also vertical half-middle 
lines connecting B % and with their neighbors in the vertical direction (B k , 
B l and B n , B m , respectively, in Fig. A. 76). As before using the periodicity 
and the divergence- free properties of u we obtain that: 



where £ v is any vertical line in Yp (Fig. A. 7b). And finally using the average 
of u and that u is zero on the inclusions we obtain: 



Proceeding similarly one can obtain that pj m — /3j n = 0. 

Hence, we obtain that the flux through any horizontal line connecting the 
neighbors equals 1, while the flux through any vertical line is zero. Also we 
note that the number of pairs for which /3y is not zero is equal to the 
total number N of the inclusions in the suspension. 

The fact that all inclusions are motionless and that there are N fluxes 
though the middle lines connecting neighbors of our network together with 
Proposition 2.1 results in a = NW = |tt^A R 5 / 2 6- 5 / 2 + 0(<r 3 / 2 ). 

B Appendix 

1 Details of construction of the trial fields 
Verification of (6.18) 

Consider the minimizer u s h of W over V s h- The test function v\ G V s u must 
satisfy the boundary conditions (6.3) on dB l and dB^ , which can be rewritten 

as 




I ui(x, y)dy for any fixed x € 





G\ = (U{ 




U{ + Ruj 1 + EtJ)g\ + R{uj 1 + <J)g\ , 



G{ = (Ul 




) 



= : (Ul 



U{ + Ruj 1 + Rj j )g{ + R(J + J)g{ . 



(B.l) 
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Hence, Vi can be written as 

v 1 = {U{ - U{ + Ruj 1 + Ruj 3 )v\ + R{uj 1 + uj 3 )v 2 , (B.2) 

where 

v\\ dBi =g[, v\\ dBi =g\, v \\ aB .=g{, v l\ gB .=gl 

Using the decomposition (B.2) we define i>i by (6.21)-(6.22). 

Next we evaluate the difference between W(vi) and W(u s h). Here denote 
£ = U[- U( + Ruj 1 + Ruj 3 and C = R{uj 1 + oj 3 ) then 

W(u sh ) = W(Cul h + Cu 2 sh ), 

where u x sh is the minimizer of W satisfying g\, g\ and u 2 h is the minimizcr 
of W satisfying g\, g 3 2 on dB % , dB 3 , respectively. 
Compute the functional W(vi) as 5 — > 0: 

W(vi) = eW(v\) + C 2 W{v\) + 2ttW(v\,v\) 

= \{U{ - U{ + Ruj 1 + R^ifn^'H^' 2 + C 2 H>?) + 2&W{v\,v{) 

=: {U{ - U{ + Ruj 1 + Ruj 3 ) 2 C\ 3 6~ 1/2 + ( 2 W(v 2 ) + 2^W(v\,v 2 1 ), 

(B.3) 

where W(v 2 ) = 0(1), W{v\,v 2 ) = 0(1) and 

C 3 = \^' 2 d^ 2 . (B.4) 

Due to parallelogram identity for any vector fields (p and i/>: 
4W(<p, ip) = W{if + t/>) - W{<p - ip), 

we obtain 

\W( Vl ) - W(u sh )\ = \W{Zv\ + (v 2 ) - W(£ul h + (u 2 sh )\ 
< ( 2 \W(v\) - W{u\ h )\ + C 2 \W(v 2 ) - W{u 2 sh )\ 

+ | \W{v\ + v\) W{u\ h + u 2 sh )\ + | \W(v\ v 2 ) - W{u\ h - u 2 sh )\ . 

(B.5) 

In order to evaluate these differences in (B.5) we construct four trial tensors 
S\ e F, k = 1, . . . , 4, such that, 

Wt(S\) < W; (cr(ul h )) , W^Sl) < W{ {cr(u 2 h )) , 

W^St) < W{ {<r{u\ h + u 2 sh )) , Wi(St) < W{ {*{u\ h u 2 sh )) , 

obtaining 
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\W{ Vl ) - W{u sh )\ < e \W{v\) - W*(S\)\ + C 2 [W{v\) - W*(St)} 

+ y \W{v\ + v\) W?(5?)| + ^ \W(v\ - v\) WZ{S\)\ . 

To construct e.g. S\ we use the following observation. If v\ is a "good" 
approximation of u\ h then we expect that <r(v\) = 2^D(v\) — p v il, where 

V 11 \\{G + F'-^G") -yG' j Vl 

with G(x) and F(x) defined in (6.22), is a "good" approximation of cr(ul h ). 
Then we correct cr(v\) so that it belongs to the set F, thus, obtaining S{. 
Now define 

ci ( G(a;)-C\ . . . 

*i - ^ G(x) _ c _ yG , (x) J m n Vt (B.6) 

with G(x) given by (6.22), and the constant C = — — -. Clearly, S] G F. 

When we evaluate the difference between W^f}) defined by (B.3) and 
W?(S\), defined by (6.13) we obtain the following estimate: 



\W(v{)-W?(S{)\=0(l), as 6^0. 



Now choose 



hence, = 0. 

We also choose <S^ = <S} and Sf — S{ and compute the corresponding 
functionals obtaining 

\W(v\ + v\) W*(S[)\ = 0(1), \W(v\ - v\) W?(S{)\ = 0(1), 

as S — > 0. Therefore, W(v\) determines the asymptotics of W(u sh ): 

iWiv^-Wius^l^E^U'iU^Lu 1 ,^), as 5^0, (B.7) 

where the error term E\ is defined in (6.18). Hence, 

W(u sh ) = [(U i -U j )-p ii +Ru i +Ru i ] 2 C^5-^ 2 +E 1 {U i ,U j ,u i ,w i ) as 6 -» 0. 

(B.8) 



Verification of (6.19) 

Consider the minimizcr u sq of W over V^ g . Choose i>2 by (6.24)-(6.25). Intro- 
duce an approximate pressure 
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Pv2 = ^{ Q J^ F( t) dt + G ' - 3 ^ 2f ') . ( B - 9 ) 



'70 

where 70 is some constant and compute 



„,„ , _ (TT i T p , ( MG' + 3y 2 F<) - p V2 n(6yF - yG" - y 3 F" 

*{v 2) - (u 2 u 2 ) ^ {6yF _ yG „ _ y3p!l) _ MG , + 3y2pl) _ pv j. 

(B.10) 

Now we construct a corrector S 2 = («S 2 )i,.j=i, 2 to (B.10) so that S 2 € F. 
Define S 2 by (6.26) where constants d (i — 1..4) are chosen so that S 2 n = 



on dFI:, , that is, 

7ij 7»j 7»j 

C 3 = -3G'( 7y ) + d^, C 4 = -9F'( 7ij ) + ^C 27 £- 

Then we compute the functional W(v 2 ) and obtain 

vk(« 2 ) - {ui - uif (cp-*/ 2 + + c^{ui - uif, (B.11) 

where 

c- = l^v 79 . c " - ^^V /a - ( B - 12 ) 



When we evaluate the difference between W(«2)> defined by (B.ll), and 
W 2 (5 2 ), defined by (6.13) we obtain the following estimate: 

\W(v 2 ) - W;{S 2 )\ < E 2 (U\ U j ), as 6 - 0, (B.13) 

where £7 2 is defined by (6.19). Therefore, W(v 2 ) determines the asymptotics 
of W(u sq ) as S — > 0: 

VK(m S9 ) = ([C/ 4 - £P] • q«) 2 (C%5- 3 ' 2 + C^- 1 / 2 ) + £ 2 ([/\ E7 J '). (B.14) 
Verification of (6.20) 

Consider the minimizer u per of W over V^, er . The test function v 3 e V^, er may 
be written in the form: 

v 3 = R(lu 1 - u j )v\ + PijRvj, (B.15) 

where 

v l\ dB i = G 3, "alas* = G i "Ham = °> = 

and the parameter flij is defined as follows. The flux of v 3 through £ij should 
satisfy / v 3 ■ nds = /?*, -Jl(Ul + U{), then 
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0H=^-^^-f vl-nds-^Ul + Ui), (B.16) 

where ^L{U[ + U{) = [ u t ■ nds. 

f S S 2 - 

Define v 3 by (6.27)-(6.28). Computing / v\ ■ nds = + we obtain 

that /3ij is denned by (B.16). 

Now we rewrite (6.27) as one vector: 

/ G(x)+3y 2 F(x) \ 
V3 -{-yG'(x)-y 3 F'(x) mIIij > (R17) 



where 

G{x) = R(lu 1 - w j )P(x) + pijRK(x), F(x) = R(uj 1 - co j )Q(x) + [3 t3 RM(x), 

(B.18) 

with P, Q, K and M defined by (6.28). 

As before to choose a trial tensor <S 3 G F per we look at the stress tensor 
cr(v 3 ) = 2[iD(v 3 ) — p V3 , where an approximate pressure p V3 can be chosen as 
in (B.9) for G, F given by (B.18). 

Then for such functions G and F define the trial tensor <S 3 <G F per by 
(6.29) where constants G\ as S2- 

Next we compute the functional W(v 3 ): 

w(v 3 ) = p% (ep- 5 / 2 + ci j s-v 2 + cis- 1 ' 2 ) 
+R 2 (oj i -uj j ) 2 ei j s- 1 / 2 

+C 5 i?V - + C 6 i?V - ut)Pij + C 7 R 2 (3 2 3 , 

where 

9 „ K/0 ,-n/2 99 „ Q /o ,_s/2 29241 -1/2.-1/2 



C y - -n,jR 5 / 2 d- 5/2 C lJ - —Tr///? 3 / 2 ^- 372 C y - 7r,,R 1/2 .r 

40 M « ' 9 ~ 2 



c« = -3^ 3/ V /2 > ei = ^R^d-;i\ a = \^r^i\ 



(B.20) 

When we evaluate the difference between W{v$), defined by (B.19)-(B.20), 
and Wg*(<S 3 ), defined by (6.14) we obtain the following estimate: 

\W(v 3 ) - W 3 *(5 3 )| < £3K,w j ,Ai), as 5 - 0, (B.21) 

where £3 is defined by (6.20). Therefore, W(v 3 ) determines the asymptotics 
of W(u per ): 
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W(u per ) = Pl (Cp- 5 ' 2 + C^- 3 / 2 + C^" 1 / 2 ) 

+R(uj 1 - uP)fc (C^'r 3 / 2 + &i8-V 2 ) (B-22) 
+R 2 (J - utyCgS- 1 ' 2 + £ 3 (wV,/%). 

2 Construction trial fields in boundary neck Ilij, (i El, j E A/i (IB) 

The construction of the trial vector and tensor fields in neck Uy connecting the 
disk B % and the quasidisk as in Fig. B.l can be done similarly to the case 
of the neck connecting two interior disks. Let B 1 be a disk near the boundary 
and £?■?' a quasidisk, neighboring to Bi. We call IP, uj % of the quasidisk B^ the 
"translational" and the "angular" velocities of the quasidisk B? , respectively. 
Now consider the neck 77y connecting the disk B z and the quasidisk B^ lying 
on the lower part of df2. The other cases are treated similarly. 




(a) 






(c) 



Fig. B.l. (a) The disk B r and quasidisk B 3 , (b) boundary conditions on B* and 
B- (c) the local coordinate system 



Recall that in general necks are not symmetric. To construct a test function 
v € Vij in the neck 27^ connecting the disk B % and quasidisk (Fig. B.l(a)) 
we extend JJ^ to 

#y = {(*> y) e K 2 : -7ij < a: < 7«, -Hij(x)/2 <y< iiy(x)/2} , (B.23) 
where 

if 7y > Ty, 

and construct w € 7J 1 (77y) such that w = C/ 1 + Rio % (n\ei — n%e.\) on <9i? 1 n 



<9/7y , v = Ax on dB^ (~l <9/7y and — / v ■ nds = /?*,- . Then the restriction 



1 

R 

of this function on iJy: v = v\n i:j would be a function from Vy (since this 
would be iJ 1 (7Ty) function and all boundary conditions would be satisfied). 
Hence, Wjj w (tt) < Wij y (t)), where it is the actual minimizer of W/^. in Vy. 

We can decompose a test function v 6 Vy (i € I, j EjVjfl B) in to the sum 
v = u t + v 1 + v 2 so that VK/Ty («) = W^. (wj.) + Wjj^ (w 2 ), and Wjjy (t> 3 ) = 
where functions Vi, v% satisfy 



58 Leonid Berlyand, Yuliya Gorb, and Alexei Novikov 

(a) u t = \{U l + U ] ) for x G dB l and x G dB^, 

(b) vi = ( U A +uj 1 (^ R2 ~ ^ ) for x G 9-8% Vi = f ^ ) lor ;r t 



•J1A 



with / v 1 -nds = p* j R- [ ut-nds^ I3* J R--[{U 1 + U r )-p 
J in hij 4 

( c ) V2 = for a; G dB\ v 2 = (^j for x G 

(B.24) 

Then the corresponding decomposition of the minimizer is: u = u t + U\ + u 2 . 

In the local coordinate system we choose these functions v\, v 2 to be as 
follows: 

U{ + yK(x) + y 2 M{x) 
K uPx-£K'(x) - ^M'{x) i 



vi=\ . 1 2 \ ' 3 V in TTj,-, (B.25) 

where 



iJ(a;) 2 if (a;) i?(x) v v ; H(x] 



2 



MW ~ H{xf H{xf + H(xy [H VU X) H{xf 



with 



Tp _ ^Ljjj _ Oil i _ o* _ ^j± T ji . „ij _ hi TT3 . n ij _ % 
AR Ul AR Ul 2 W -/ *« 4/^ P 4i?^ P 2 



and 



where 



ax + yF(x) + y 2 G(x) 

V2= {ui-a y -y-lF'(x)-y-l G '(x)) in77 - (B ' 2(,) 



8ax: 24(t/j — l/|)ar 12ax: 48(^ - 

[X) ~~Wx) i^F ' {X) ~~H~W + H{xf ' 

We construct a lower bound for the viscous dissipation rate in the neck 
near the boundary we consider a dual variational problem in a symmetric neck 
(B.23) with 



Hi 



,7ij: if 7y < 7ij; 



that is, the problem of the maximization of the dual functional (6.13) with G\ 
and G 3 given by (B.24a) in the dual space Fper given by (6.16), all defined 
over the neck 77^ . Indeed, 

W nij {u)>W h (u)>Wl{S). 
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Here the test tensor S is decomposed into the sum S = Si + S 2 , with Si G 
F P er, S 2 G F corresponding to two trial vector fields V\, v 2 defined above. 
We choose 



Si 



-2 f* o M(t)dt + Cix K + 2yM - Ciy - C 2 
K + 2yM - Ciy - C 2 -yK' - y 2 M' 



in 77, 



(B.27) 



where constants C\ — 2M(jij), C 2 — K(jij). This tensor is divergence free 
in Ilij , and on the lateral boundary dll^ of the neck 77^ it takes values: 

X + - -2 f % ' M(x)dx + Cijij, X + = -2 I 7 " M(t)dt - Cijij, 

J la J la 

hence Si G F per . 

Now we compute the functional W(vi) in the neck connecting a disk and 
a quasidisk: 



W(vi) = pi 



-3/2 



20889 



7rnR 1/2 5-/ /2 



- 2,TT^R 1/2 l3 i0 Ruj l 5~ 112 + 4tt^ 1 / 2 (u[ - U{ + Rlu'Y 5^ 1/2 

- ^iiR 1/2 (u{ - U{ + Ru^ R(lo 1 - u j )6Tj 1/2 

+ K^ 2 R 2 (lu* -<J)% 1/2 + E b 1 (U i ,w i ,f,0 ij ) 



3/2 



(B.28) 



= : pl \e[ j 6- 5 / 2 + ep- 3 / 2 + cp-v 

+ difci (U{ - U{ + Raj 1 ) (T 3 / 2 + CiP^Riuj 1 - u j )6 

+ C£% (ill - U{ + Rw^ 8- 1 ' 2 + Cyp i3 R{u l - u^S- 1 ' 2 

+ &i /JyitoVT 1 / 2 + Cl j (U{ - U{ + Rw i y 5- 1 ' 2 

+ Cio (yi - u ( + - u J )5- x / 2 

+ C\\R 2 (^ - uPfS- 1 ' 2 + E b i(U\ u\ f, fa), 
where the error term 

E\{U\uj\ f,fa) = /x (C 7 /3 2 - + C & [(U l - U j )-p^} 2 + C 9 R 2 (^) 2 + C 10 i?V) 2 ) 

(B.29) 

with universal constants Ck, k = 7, ... ,10 and 
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C« = 67 r^ 3/2 rf- 3/2 , C« = -torfl*d£'\ Ci = f^R^\ 

eg, = -touBWdJ", = 

When we evaluate the difference between defined by (B.28), and 

W 3 *(Si), defined by (6.14) with the functions G 4 fe , G?£ given by (B.24a) we 
obtain the following estimate: 

\W{v x )-W^S x )\ <£j(t/W, /,/%), as 5^0. (B.30) 

Therefore, W(«i) determines the asymptotics of W{u\) in the neck connecting 
the disk and quasidisk: 

W( Ul ) = 0* \&l5- b ' 2 + Bp- 3 / 2 + Bp- 1 / 2 

+ BiiPij [(C/ 4 - W) ■ p ij + Rcj 1 ] 6~ 3 / 2 + BgPaRiw* - ^)S- 3/2 

+ [{U l - U j ) ■ p ij + Rlj 1 ] 5~ 1 ' 2 + ByfaR^ - uj j )8-^ 2 

+ B'iH^Ruj^- 1 ! 2 + Bl j [{IP - U j ) ■ p l > + Ruj 1 ] 2 5- 1 ' 2 

+ B% [(IT - U j ) ■ p l] + Ruj 1 ] R(uj 1 - (J)S~ 1/2 

+ B\\R 2 (lu* -uifS- 1 ' 2 + E b 1 (U i ,w i ,f,/3 ij ). 

(B.31) 

Now we choose S 2 € F as follows: 

(«S 2 )n = fj, ^3yF' + 3y 2 G' -2 J G + y (C 2 + C 3 y + C A y 2 ) 
-(C 5 + C 6 y + C 7 y 2 ))'! 

(S 2 )i2 = (5 2 ) 2 i =fi(F + 2yG- ^-y 2 F" - y 3 G" - C lX - C 2 xy - C 3 x^- - C 4 x 



(S 2 ) 22 = V {-yF' y 2 G' - 2 J" G + \y 3 F<" + \y 4 G" 



'7o 

y 2 y 3 y 4 

+Ciy + C 2 ^ + C 3 y - + C i ^ 

(B.32) 

where constants d (i = 1..7) are chosen so that S 2 n = on dFlf-, that is, 
C^—F^), C 2 = ^-G( lij ), C 3 = -^-F"( lij ), C t = -— G"( 7ij ), 
C 5 = C^-f - 2 jT" G, C 6 = C^-f + 3F'( 7y ), C 7 - cjf + 3G'( 7ij ). 
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We now compute the functional W(v 2 ): 
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w(v 2 ) = {ui - uif 

+ GnfiR 1 / 2 ^ - Ui)Ra5; 3 1/2 + E h 2 {U\ f) 

+ e[{(Ul - U^RaS- 1 ' 2 + E b 2 (U\ /). 
where the error term 

E b 2 (U i ,f)=»C 11 [(U i -U j )-q i i} 2 , 
with universal constant Cn and 



Cg = Z^R^d^' 2 , C% = ^Ri/ 2 d^ 2 , Cg = 6^rV% 



(B.33) 



(B.34) 



1/2 



When we evaluate the difference between W(v 2 ) 1 defined by (B.33), and 
W 2 (S 2 ), defined by (6.13) we obtain the following estimate: 



\W(v 2 ) - W;(S 2 )\ < E b (U\ /), as S - 0. 



(B.35) 



Therefore, W(v 2 ) determines the asymptotics of W(u 2 ) in the neck connecting 
the disk and quasidisk: 



W{u 2 ) = [(IT - U j ) ■ q^] 2 (Bgs- 3 / 2 + sgs- 1 ' 2 ) 

+ B\{ [(IP - U j ) ■ q ij ] RaS- 1 ' 2 + E b (U\f). 



(B.36) 



C Appendix 

1 Positive definiteness of Qij in (2.36) 

Lemma C.l. The quadratic form Qij, defined by (2.36), (A.l), is positive 
definite. 

We show the positive definiteness of Qij with coefficients (A.l) for i g I, 
j € Afi n I (the case when j ejV.nl can be treated similarly) . To this end 
we just need to show that Qij(U l , C/ 3 , lo^ ,(3ij) can be written as a sum of 
squares. Indeed, 
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207 
320 

. , , 249 , 13491 \ 
+ %{ 2S + W + 17920^) 
T3_. , 

+ 

+ E ij (U\U>,u;\w>,f3 ij ) 



-R(J - uP) - ( -6- 1 - " 
4 v ' y 3 \ 2 16i? 



^0= 1 1 [(17* - U*) ■ p« + iiu/ + ft^] 

+ [(C/*-^).,«] a (|iM- 1 + | 

„, / fir , 249 \ 2 645273 . 

>^-')-"«(k'-Ti)] 2 } 



+ E ij {U\U\^,w\(3 ij ) 

2 Existence and uniqueness of the solution to (2.4) 
Lemma C.2. There exists a unique minimizer ofWn F {-) over V. 

Proof. First of all we note that V is not empty. Indeed, let w be a solution 
to Stokes equation satisfying the Dirichlet boundary conditions, that is, of 

' (a) /iAtt = Vp, x G ftp 

(b) V • u = 0, xe n F 
... (CI) 

(c) u = J7q + RLo l {n\e 2 - n\e{), x e dB\ i=l,...,N 

K (d) u = f, x € dfl 

with some given U l Q and wj, i = 1, . . ., N, which is known to exist. Then 

w e V. 

The existence of the minimizer u would be shown in few steps. 

Step 1. Set m := inf Wq f (-). If m = +oo then we arc done. So we hence- 
forth assume that m is finite. Select a minimizing sequence {uk} ( ^L 1 of the 
functional Wq f {-), that is, Wn F (uk) — > m as fc — > oo. 

5"tep Choose a sequence {uk}'^L 1 , where each xik solves a Stokes equation 
(C.l) with boundary conditions (C.lc) on dB % which are the same as the 
corresponding boundary conditions of u^. Hence, 

m < Wn F (u k ) < Wq f (u k ) — > m as k -> oo. 
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Now we claim that there exists a convergent subsequence °f values 

of velocities on each disk, that is, there exists ip l such that ^sL — > tp 1 as 
j — > oo on each dB % , i = 1, . . . , N 

We prove this assertion by contradiction. Assume there exists such a disk 
B ia and a subsequence {<pl)} = {Ufy + 17^} on it that y>|°. -> oo as j -> 
oo, where U 1 ^ is the translational velocity and fl^, is the rotational part. 
We select a subsequence {ukj} £ V of the above sequence that is, 

{ukj} solves the Stokes problem with fixed values ip^ on the disks B % . Since 
Wq f (Ukj ) is bounded and cp]l'j — > oo we have 

where e]y is a small arbitrary number. 

Without loss of generality, we assume that the disk B l ° lies in the corner 
of the domain Q : that is, near two parts of the boundary dft (see Fig. C.l), 
where values of Ukj are prescribed and equal to /. 



an 




Fig. C.l. 



First consider a neck iTj j' and choose a symmetric neck n ioi i containing 
TT^j^, that is, n ioi i o D n ioi ' Q . For simplicity we drop indices fcj indicating the 
subsequence and zq indicating the disk number. Then denote ip = (ip±,ip2) — 



(U, O) = (t/i +n u U 2 + n 2 ) and 



(it, v). Choose a point (xq, yo) on the 



lower part of dfl and (x,y) on dB l " (Fig. C.l). By integrating along the line 
connecting x and xq we obtain: 



v(x,y) = 



f2(x ,y ) 



m 



v dv 



or 



|w(x,y)| < r/ 2 + 



yo d V 

"(7 i0i& )/2 ^ 

dy 



(x,y)dy, 



dy, 
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where v(x,y) — <p 2 and r\ 2 is some small number since ^%^p^ — * as j ' — * 00 

and H("fi a i> )/2 is a half-hight of the neck n igi > g . Then integrate over the arc 
a^ o of <9-B 10 which is an upper boundary boundary of the neck 2Jj j' : 



/ uds 




Ja 1 

'0 





\v\ds < \al \r) 2 



n 



9w 
dy 



dx 



and 



■jds 



< Ci772+C 2 



0*0 



dv 

dy 



(x,y) 



dx < Cir)2+C 2 Wn F {u) < Cn] 2 +C 2 e. 




Consider 

/ vds= [ (U 2 + n 2 )ds = \a\\U 2 + ( Q 2 ds. 

Ja 1 Ja 1 Ja 1 

ZQ % l Q 

Note that J 1 Q 2 ds = since the second components of the vector fi cancel 

each other when we write an integral as a sum over infinitesimal elements of 
the symmetric arc a\ o (see Fig. C.2). Hence, 

CW 2 2 < Cir] 2 + C 2 e, 

where r/ 2 = (^2)^ and e = e]y are small arbitrary numbers. So the second 
component of the vector IT£. tends to zero as j — > 00. 

Similarly, one can show that the first component of the vector U]y tends 
to zero as j — > 00 considering the neck n iQi ». 
Thus, we prove that U 1 ^ — > as j — > 00. 

To show that fl^ — > as j — > 00 we consider only the left-hand side of 
the neck 11^ . Repeating all computations described above, we obtain that 
the integral of D 2 over the left half of the symmetric arc a^f tends to zero as 
j -> 00: 
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/ Q 2 ds < Ci77 2 + C 2 e, 

Ja ie 

'0 

But l?2 = cos 9, where u>^- is a constant and 9 is an angle of the vector 
fi 1 ^ with vertical axis. Thus, 9 > everywhere in that half of the neck, hence 
J? 2 ~~ * 0. Similarly, one can show that Q\ — > 0, hence, the whole subsequence 
(p^ is convergent what contradicts our assumption. 




Fig. C.3. 



Note that if B ta is not near dfl one can use the same argument using two 
nearest neighboring disks B rni and B m2 of B l ° where there exist convergent 
subsequences if]™ 1 and <p™f, respectively. We can always assume that since 
our graph is connected, that is there exist paths connecting the given disk 
with all parts of the boundary dQ. 

Step 3. Denote the limit of convergent subsequences on each disk B l by 
v l , i — 1, . . . , N and consider the subsequence {ukj} G V such that on the 
boundary of each disk dB % : Ukj — > v l as j — ► oo. 

Step 4- Choose a sequence {uj} £ V such that each Uj solves the Stokes 
problem (4.5) with iij = v l on dB % , i = 1, . . . , N. 

Step 5. Note that 

Wn F (uj) < Wn F (u kj ) + W S2f (%) + 2^W^^W^~), (C.2) 
where Uj — u^j ~ Uj , that is solves 

' (a) fiAu — Vp, x £ flp 

(b) V • u = 0, a; £ Q F 

< 

(c) m = v l kj - v\ x £ dB\ i = 1, . . . , N 
k (d) w = 0, x£dQ 

Hence, (see e.g. [26]) ||%|| H i (f2F) < C||«^. - f l || ff i/2 (asi) -> as j -> oo. 
Then from (C.2) we obtain 

Wftj^tlj) — > m as j — > oo. 
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Step 6. Choose u = Ui and u°- = Uj — u where vP- € Hq{Q f ) such that 
Wq f {uj -\- u) — > in as j ► oo . 
Step 7. Consider 




Here we used the integration by parts (all boundary terms disappear since 
w° G Hl{Qp) and the divergence free property of Thus, continuing (C.3) 
we obtain 



Ci\\Wu°\\ L , { n F ) +C 2 < C 3 Wn F (u° 3 ) + C2 

So the sequence iij is bounded in L 2 (f2p), as well as 

HVSj-IU^) < ||V Mj °|| L2(r2F) + ||V5|| i2(flp ) < C 5 . 

Note that u € H 1 (J7_f) as a solution to the Stokes problem [26]. 

Therefore, from all the above Uj is bounded sequence in H 1 (f2p). 

Step 8. From the previous step it follows that there exists a vector field 
u e H 1 (Hp) such that Uj — 1 u (weakly) in H 1 (Hp) as j — ► 00. 

We have to show that u € V. Recall that for u e V: iij — u e Hq{Qp), 
which is closed linear subspace of H 1 (f2p), hence, Hq(Qp) is weakly closed 
(Mazur's theorem), thus 

iij — u ^ u — u, 

therefore, u e V. 
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Step 9. Finally, we show that Wq f {u) — vn. The functional Wq f {-) is 
convex in Vii hence it is weakly lower semi-continues, that is, 

Wn F {u) < liminf Uj = m. 

Since u e V, we have 

W( 2f {u) =m = imnWn F (-)- 

3 Existence and uniqueness of the solution to (2.10) 
Lemma C.3. There exists a unique minimizer ofWn(-) over Vn- 

Proof. To show the existence of a minimizer of functional Wn(v) for v e Vn 
we have to prove that it is coercive and convex in the gradient of v and that 
Vn is nonempty convex and closed set. 

Step 1. Set m := infy^ Wn- If m = +oo, we are done. So we hereafter 
assume that m is finite. Select a minimizing sequence {u n }^i C Vn of 
Wn(-), that is, Wn(u n ) — *• m as n — > oo. 

Step 2. By the argument similar to Step 2 of Lemma C.2 one can choose 
a bounded subsequence of the velocities (which are constant vectors) on the 
disks B l , i — 1, . . . , N. Denote such a subsequence by {v l n }^ =l7 i — 1, . . . , N, 
that is, there exists a positive Co such that for all n and every i = 1, . . . , N: 
\v l n \ < Co- Denote the limits of those subsequences by v l . 

Note that by the Korn's inequality (see e.g. [18]) we have: 

||Vtx n ||| 2(J7) < C n (W n {u n ) + \\u n \\ 2 mn ^ , (C.4) 

where Cn is some constant. Therefore, we first show that the sequence u n is 
bounded in L 2 (II). 

Step 3. Select a sequence w n of solutions to the Stokes problem (C.l) with 
w n = v\ on each dB l , and the functions w that solves (C.l) with w — v l on 
dB\ i = l,...,N. Then, 

||w n ||ffi(fi F ) < \\w n - w\\ H i(n F ) + \\w\\ h i {s - 2f ), 

where both norms of the right-hand side are bounded [26] as solutions to 
Dirichlet problem, thus, ||i« ra ||i? 1 (J7) < C. 

Define a sequence tp n = w n \n — u n in iT, which is zero vectors on each 
sphere dB l and zero vector on the external boundary dfl. Using (C.4) for the 
sequence ip n and the fact that 

Wn(<fi n ) < Wn(w n ) + W n (u n ) + 2^W n {w n )W n (u n ) 
<C + (m + e„) + 2v/C(m + e„) < Ci, 

we prove that (p n is bounded in L 2 (II). 
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Consider the neck 77^- (Fig. C.4) between adjacent disks B l and BK 
Choose a point (x ,yo) on the sphere dB^ and (x, y) € 77^-, and consider 

<P 1 n(*,V) = J'^&VW' 

for almost all x, where the superscript here indicates a component of the the 
vector (p n = (y>* , ip 2 ^). Then we obtain the following estimate: 



\<fl{x,y)\ 2 <C 2 [ 

Jx 



Hh)/2\Q 1 2 



and 7 = max(7 i -,7^). Now we integrate the last inequality with respect to y 
between 7^ and 7^, and with respect to x along the arc a, obtaining 



I 



^ n (x,y)\ 2 dx = WvUx^Wl^ < C 2 \ 7 ± -7«l#(7)/2 = C 3 . 

o 1 

Here we recall that J n \-^\ dx is bounded from above by Wn(<P n )- So the 
first components of <p„ form a bounded sequence in L 2 (77y) for any neck 

/7 tf g n. 

To show that the second components of ip n is a bounded sequence we 
consider another neck 77^ (see Fig. C.4(b)) that forms a small angle a with 
the direction of the second component <p 2 n . Then by triangle inequality we 
have 

||¥>nlU 2 (iI ifc ) < C ( a ) (HVnlU 2 (iT 4fc ) + ||<PnlU 2 (iT ifc )) i 

where the constant C(a) depends on the small angle a. Since the viscous 
dissipation rate is invariant under rotation we evaluate the norm of ip\ to 
obtain an estimate for (p 2 n . 

In the system of coordinate with y-axis along the neck 77^ we consider 



<fl{x,y) = J ^(x,r])dr], 



where the point (x ,yo) lies on dB k . As before, integrating with respect 
to x over (7^,7^) and with respect to y over the arc at and using that 

In \lfif\ 2 d x i s bounded from above by Wn(f n ), we obtain the following 
estimate: 



<pl(x,y)\ 2 dx = yl(x,y)\\ 2 LHnii) < C 4 . 



Hence, the sequence <p„ is bounded in L 2 (II). In the view of the Korn's in- 
equality (C.4) we obtain that (p n is bounded in 7? 1 (JT). By triangle inequality 
we obtain that sequence u n is bounded in 77^(77) as well. 

The boundedness of u n in 7J 1 (JT) implies there exists a subsequence u nn i 
and a function u e 77 1 (77 i:) ) such that u nn i — u (weakly) in 7J 1 (77 i:) ) as 
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Step 4- Vn is convex if for any it, v e Vn their convex combination is 
from Vij : w = ( 1 — X)u + Xv e V/j, < A < 1. To show that w; €E Vn" we have 
to prove V ■ w = in II, w = f on 977 n <9i?, u> = IP + Rw l (n\e2 — n 2 ei) 



on and / «; • nds — 0, which all clearly hold. 

JdA iih 

Step 5. To prove that Vn is closed we consider a sequence {v n }%L 1 , v n S 
V/j for all n and u n — > v in iJ 1 (J7), and show that t> S V/j as well. For 
this we need to check that V ■ v — in II in the weak sense, and v = 

U' 1 + RLo l (n[e 2 - n\e x ) on <9_B', i = I, . . . , N, / v • nds = 0, and w = / if 

JdA iJk 

dllndil, which understood in the sense of trace. Any convergent sequence in 

if 1 (77) is weakly convergent, hence, / v n ■ Vtpdx — > / v ■ Vtpdx for every 

Jn Jn 
(p G Cq°(II). Thus the incompressibility of v holds. 

The boundary conditions of the definition of Vn are considered in the 

sense of trace, that is, for the H 1 function v n there exist a linear bounded map 

7(v n ) and a constant C such that ||7(vn)||L 2 (aii) < C||u n ||iji(jj). Therefore, 

we have 

^2 / (70fc) - l(v)) ■ nds = ^2 / j(v n -v)-nds 

< \ dA ijk\ 1/2 h(v n - «)lli / 2 2 (9Ai . fc) < CihK - «)lli 2 2 (ai7 ) 

1/2 

< Call 
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t V /2 

where C 2 = \C \ dA ijk\ , and dll = \JdB l U (J dA ijk 

\ iei,j,keAfi J iei jei,i,fceA/'« 

The remaining boundary conditions can be proved similarly. 

Finally, we show that Wrr(-) is convex in Vu. For this we write 

Wjt(Vm) = n / L(Vu)dx, 
Jn 

where L(p) = p 2 n + \{p\2 + P21) 2 + V22 and P = (Pij)i,j=i,2 is 2 x 2 matrix. 
Recall that by definition convexity means that 

2 

X! L Pi 3 Pki(p) a ij a ki > 0, for any A = (a y ) 4j=1 , 2 . 

i,j,k,l—l 

Therefore, we have 



L PllPll 


= ^P22P22 


= 2, L pl2Pl2 — L P2lP21 — L P12P21 — L P2lPl2 = 1, 




^P11P21 


— L piip22 Lp 12pil — L pi2 p 22 Lp 2ipil — Lp. 2l p 22 




= ^P22P12 


= Lp 22P2l = 0. 


Hence, 







2a\ l + 2a\ 2 + a\ 2 + a% l + ai 2 a 12 + a 2 ia 2 i = 2af 1 + 2a% 2 + (012 + fl2i) 2 > 0, 

thus, Wn(-) is convex in Vu. This implies that Wn is convex, it is weakly 
lower semicontinuous: 

Wn{u) < liminf Wn{u n ) = m, 

n 

and Wn{u) — m. 

By Mazur's Theorem Vn is weakly closed, hence, u G Vn- 
Step 6. To show a uniqueness of the minimizcr we proceed by contradiction. 
Assume that there exist two minimizcrs of this problem u and u. Consider 
their difference w = u — v that solves (4.5) with / = on d£2 and some 

W % + Rijj]^{n\e2 — n 2 ei) on dB l , i = 1, . . . , N. By multiplying equation (4.5a) 
by w and integrating by parts over II one obtains that Wn{w) = 0. 

Using the argument given by Step 2 of Lemma C.2 we can show that w = 
on each disk B l , i — 1, . . . , N and the argument of Step 3 of the current lemma 
we can prove that w — in II. 

4 Existence and uniqueness of the solution to (4.4) 

Suppose there are two solutions to (4.4). Then their difference satisfies 
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(a) fj.Au = Vp 

(b) V • u = 0, 
(c') « = 0, 
(c") u = 0, 



3? g n^j , 



w • nds = 0, 



(C.5) 



(e) er(w)n = -pfjU, x e dll*, 
(/) u = 0, x e 577^ n <9J?. 

Integrating by parts the incompressibility condition (C.5)(b) over half-necks 
and using (C.5)(b),(c'),(c"),(d),(e), we obtain: 



/ u ■ nds = I 

Jan+ Jai 



u ■ nds = 0. 



Multiplying (C.5) (a) by u and integrating by parts, we obtain 

Wnij (u) = I u- cr(u)nds = p + / u ■ nds + p~ u ■ nds = 0. 
Jan zi Jan+ Jdnr. 

Hence u must represent a rigid-body rotation. Due to our homogeneous 
boundary conditions, this motion must be identically zero. 

Note, that since u the solution to (C.5) must be identically zero, it implies 
that pfj — 0. Hence we showed uniqueness of pf- as well. 



5 Dual problem to (2.10) 

Lemma C.4. Let u per = {u\ er \u\ er ) is the minimizer ofWn i:j {-) over V per . 
Then the stress tensor er(u per ) is the maximizer ofW^(-), defined by (6.14), 
over the set F per , defined by (6.16). 

Proof. For simplicity we use a notation <x for <r(u per ). To prove that <x 
is a maximizer of WS (•) we consider an arbitrary tensor T £ F per , such that 
Tn = XT n on wnere Xt' Xt are arbitrary numbers, then evaluate 

W£(cr + T) and show that W 3 *(cr) > VK 3 *(cr + T). 



Denote the first linear part of (6.14) by 



£(<r) = 



H (%) Jan+ 



n ■ crnds 



n ■ crnds 



an: 



n-ij nds 



R(3*j 



1 



1 



R0ij(Pij -p«). 



while 
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C(T) ^ 



n-Tnds , 



n ■ Tnds 



Also denote 

M{(t) = R^-uP) 



G\ ■ crnds + 



as* 



/ 

JOB? 



Gg • crnds 



Now we consider the quadratic term of (6.14) and notice that for er + T 
one has 



4m Jm 



(cr + T):(cr + T) 



tr(<r + T) 2 



1 



cr : cr + 2er T + T T 
= w(cr) + 2w(cr,T)+w(T), 
where 



(tax) 



dx 

2 



trertrT 



(tvTf 



W((T,T) = - 



1 



^ in,, 

and the corresponding quadratic form: 



(cr:T)- 



tr er tr T 



dx, 



dx 



(C.6) 



w( 



(er:*) 



(trrrf 



dx. 



Hence, W£ (cr + T) can be written as: 

W; (cr + T) = C(cr) + C(T) +M(cr) + M(T) + w(a) + 2w(er, T) + w(T), 

where W£(er) = C(cr) + A4(cr) + w(cr) in these notations. 
Taking into account that 



cr = cr(u per ) 



dy + dx ) Z A* dy Vo 



and 



we consider 



Pa = -- 



trer 

2~' 
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trertrT" 



(<r:T) 



dx 



1 

2 M im 



2/i 



a< 



Pa Tn + 2/Z 



du\ er du\ 



dy dx 



Ti 2 



d 



per 



dx 



2 V Jn zj 



r\ per o per r\ per o per 

2»^—T n + ^^—T 12 + 2/i— ^T 12 + 2/i— ^T 22 
ax ay ax ay 



dx 



n per o per n per <~i per 

du{ du p du p 2 0u p 2 

■in H ^ ii5 H ^ -112 i ^ I; 



ax 



ay 



dx 



dy 



22 



(fa 



p er 0T\\ p er 0T\2 per ^-^~12 . per ^-^~22 



dy ' 1 dx ' 2 ay 
- / (uf er niTn + u per n 2 Ti2 + iF 2 r n\T\ 2 + u P f r n 2 T 2 2)ds 



enz 



■Tnds 



Hiipey ' 'J- * t\i ds 



-R(lu 1 - uj j ) (I G\ ■ Tnds + [ Gi- Tnds] 



(C.7) 



where the volume integral over 77^ disappears due to divergence-free property 
of the tensor T. Continuing (C.7) and taking into account Tn — Xr n on 

and the fact that u per satisfies f u per ■ nds = ±/3*j (this is from the 



ianf 



divergence free condition of u per in 77*) we obtain: 

2w(cr, T) = — Xt I u per ■ nds — Xt / u per -nds 
Jan+ Jan- 

-R(uj 1 -u j )( [ Gl- Tnds + f Gi. ■ Tnds] 

KJdB* JdBi J 

= -X$Rfti + XtRKo - M(T) = -C(T) - M(T). 
Finally, we observe that 



w(T) = - 



1 

4/t Jn 



(T : T) 



(trT) 2 



dx = — - 



1 

Jn t 



1 



r(T 11 -r 22 ) 2 + 2T 1 2 2 



dx < 0. 



Therefore, from all the above we conclude 

W^(cr + T) = C(a)+C(T)+w(cr) + 2w(a,T)+w(T) 

= W 3 *(<t) + C(T) - M{T) +M(T)- C(T) + w{T) 

= W;{<t)+w{T)<W;{ct), 

and this inequality holds for any T <G F per . Thus, cr{u per ) is a maximizer of 
the functional (6.14). 
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